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COMPUTATIONS  OF  NONLINEAR  GRAVITY  WAVES 
BY  A  DESINGULARJZED  BOUNDARY  INTEGRAL 

METHOD 


by 

Yusong  Cao 

/ 


Chairpersons:  Robert  F.  Beck,  William  W.  Schultz 

A  desingularized  boundary  integral  equation  method  combined  with  an  Eulerian- 
Lagrangian  time-stepping  technique  is  developed  for  nonlinear  gravity  wave  prob¬ 
lems.  The  desingularization  distance  between  the  boundary  and  the  sources  is  related 
to  the  local  mesh  size  to  ensure  convergence.  Tests  for  some  simple  problems  show 
that  desingularization  significantly  reduces  the  computer  time  required  to  compute 
the  influence  matrix  of  the  resulting  algebraic  system.  The  algebraic  system  is  still 
adequately  well-conditioned  to  allow  fast  iterative  solutions.  Accurate  solutions  can 
be  obtained  for  a  large  range  of  desingularization  distzuices  on  the  order  of  the  mesh 
size. 

Several  nonlinear  water  wave  problems  Me  then  investigated.  The  first  prob¬ 
lem  considers  upstream  runaway  solitons  due  to  a  disturbance  moving  near  critical 
speed  in  two-dimensional  shallow  water.  Results  from  the  desingularized  method 
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with  the  fully  nonlinear  free  surface  boundary  condition  agree  well  to  those  using 
the  fKdV  model  for  weak  disturbances.  The  fully  nonlinear  model  predicts  larger 
solitons  than  the  fKdV  model  for  strong  disturbances  and  also  predicts  the  breaking 
of  waves  for  some  stronger  disturbances.  Next,  the  problem  of  three-dimensional 
waves  due  to  a  submerged  moving  spheroid  show  good  comparison  to  those  from 
other  algorithms.  Finally,  the  generation  of  inner-angle  wavepackets  in  the  wake  of 
a  ship  is  investigated.  The  three  most  probable  causes  of  the  wavepackets  are  exam¬ 
ined:  interference  of  the  wave  systems  by  the  bow  and  stem;  free-surface  nonlinear 
effects;  and  wake  unsteadiness  due  to  translation  and  os'~.illation  of  the  disturbance. 
The  wake  is  studied  with  nonlinear  calculations  using  the  desingularized  method 
and  with  linear  calculations  using  a  time-domain  Green  fimction  and  the  stationary 
phase  method.  It  is  shown  that  nonlinear  effects  are  not  essentisJ  to  the  generation 
and  persistence  of  inner-angle  wavepackets;  the  phenomenon  can  be  explained  by 
unsteady  linear  theory. 
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CHAPTER  I 


INTRODUCTION 


Problems  involving  nonlinearity  of  free  surface  waves  in  hydrodynamics  have  be¬ 
come  of  increasing  interest  since  many  wave  phenomena  cannot  be  explained  by  lineair 
theories.  For  instance,  in  the  generation  of  the  solitons  running  away  upstream  of  a 
disturbance  moving  near  the  critical  speed  in  shallow  water,  the  nonlinear  effect  plays 
an  essential  role  (Wu  k  Wu  1982,  Ertekin,  Webster  k  Wehausen  1986).  Another 
phenomenon  is  the  observation  that  a  ship  travelling  in  deep  water  may  generate 
a  relatively  steep,  isolated  wavepacket  within  the  diverging  portion  of  the  Kelvin 
wake  and  the  packet  persists  for  a  very  large  distance  behind  a  ship.  Although  the 
mechanism  responsible  for  this  wave  feature  is  not  entirely  cle£tr,  it  is  believed  that 
nonlinear  effects  Me  important  in  the  generation  and  persistence  of  the  so-called 
inner-angle  solitary  wavepacket  (Brown  et  al.  1989,  Akylas,  Kung  k  Hall  1988  and 
Hall  k  Buchsbaum  1990).  Many  aspects  of  other  wave  problems,  such  as  the  wave 
resistance  of  a  ship  and  the  motions  of  a  ship  in  waves,  require  accurate  computations 
accounting  for  nonlinear  effects. 

Since  there  are  very  limited  problems  for  which  analytic  solutions  can  be  easily 
found,  a  numerical  method  is  necessary.  Two  major  difficulties  encountered  in  the 
study  of  nonlinear  waves  are  that  the  free  surface  boundary  conditions  are  nonlinear 
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and  the  free  surface  boundary  is  unknown  o  priori.  Perturbation  methods  are  most 
widely  used  to  consider  nonlinearity  by  including  higher  order  terms  in  the  free  sur¬ 
face  bound2Lry  condition.  However,  the  boundary  condition  is  still  satisfied  on  the 
undisturbed  flat  free  surface.  Another  alternative  is  to  use  an  iterative  procedure. 
For  example,  in  the  wave  resistance  calculations  by  Jensen,  Mi  &  Soding  (1986),  the 
boundary  condition  is  satisfied  on  a  deformed  free  surface  using  an  iterative  proce¬ 
dure.  In  each  iteration,  a  lineairized  bovmdary  condition  based  on  the  solution  from 
the  last  iteration  is  satisfied  on  the  free  surface  determined  from  the  last  iteration. 
Their  method  is  not  suitable  for  unsteady  waves. 

The  use  of  fully  nonlinear  and  time-dependent  free  surface  boundtiry  conditions 
in  numerical  computations  was  introduced  by  Longuet-Higgins  &  Cokelet  (1976) 
using  a  mixed  Eulerian-Lagrangian  method  for  two-dimensional  surface  waves.  This 
procedure  requires  at  each  time  step:  1)  solving  a  boundary  value  problem,  which 
itself  is  linear,  in  an  Eulerian  frame  to  obtain  the  velocity  on  the  free  surface;  and  2) 
updating  the  free  surface  points,  which  construct  the  free  surface,  by  integrating  the 
nonlinear  kinematic  and  dynamic  free  surface  boundary  conditions  with  respect  to 
time.  Variations  of  this  method  have  been  used  for  a  variety  of  nonlinear  free  surface 
problems  in  two  dimensions  (Vinje  &  Brevig  1981  and  Baker  1982)  The  method  can 
be  easily  extended  to  three-dimensional  problems.  For  exaunple,  the  method  was 
used  by  Dommermuth  &  Yue  (1987)  and  Kang  (1988)  to  solve  several  axisymmetric 
nonlinear  free  surface  wave  problems.  Cao,  Schultz  &  Beck  (1990,  1991a,  1991b) 
used  this  time-stepping  procedure  combined  with  a  desingularized  boundary  integral 
technique  to  study  nonlinear  waves  caused  by  disturbances  moving  tmder  the  free 
surface. 
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The  discretization  of  the  free  surface  is  necessary  in  numerical  calculations.  For 
the  time-stepping  procedure  to  be  of  practical  use,  it  is  important  to  have  an  ef¬ 
fective  solution  method  for  the  botindaxy  value  problem  since  it  requires  most  of 
the  computation  time.  A  boundairy  integral  method  is  powerful  because  it  reduces 
the  computational  domain  by  one  dimension.  In  conventional  boundary  integral 
formulations,  the  singularity  of  the  fundamental  solution  is  placed  on  the  domain 
boundary.  This  requires  special  evaluation  of  singular  integrands,  which  can  result 
in  costly  numerical  calculations. 

If  the  singularity  of  the  fundamental  solution  is  placed  away  from  the  boundary 
and  outside  the  domain  of  the  problem,  a  desingularized  boundary  integral  equation 
is  obtained  (von  Karman  1930,  Kupradze  1967,  Webster  1975,  Heise  1978,  Schultz 
&  Hong  1989  and  Cao  et  al.  1990,  1991a,  1991b).  One  immediate  advantage  of 
the  desingularization  is  that  the  computational  effort  required  for  the  evaluation  of 
the  integrals  can  be  greatly  reduced.  However,  there  is  a  question  about  how  to 
do  the  desingularization,  or  more  precisely,  how  to  select  a  proper  desingularization 
distance. 

In  this  thesis,  the  desingularized  boundary  integral  method  is  re-examined  and 
some  of  the  above  mentioned  nonlinear  wave  problems  aire  investigated  using  the 
desingularized  method.  A  brief  synopsis  of  each  of  the  following  chapters  follows. 

In  Chapter  2,  the  initial-boundary  value  problem  is  formed  for  nonlinear  free 
surface  wave  problems  under  the  irrotational  flow  assumption.  The  solution  proce¬ 
dure  for  the  problem  and  the  numerical  implementation  of  the  method  are  briefly 
described  for  waves  generated  by  a  free  surface  pressure  disturbance  or  underwater 
disturbances  {t.g.  simple  singularities  and  a  submerged  body). 
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Chapter  3,  representing  an  incorporation  of  Cao,  Schultz  &  Beck  (1991a),  re¬ 
examines  the  desingularized  boundary  integral  method  for  some  simple  potential 
problems  with  the  boundaries  extending  to  infinity.  The  uniqueness  of  the  solution 
and  the  convergence  of  the  numerical  solution  to  the  desigulzurized  integral  equation 
are  discussed.  A  formula  to  determine  the  desingularization  distance  is  proposed. 
The  effects  of  the  desingularization  distance,  the  mesh  size  and  the  truncation  of 
boundary  are  tested.  The  method  combined  with  Eulerian-Lagrangiam  time-stepping 
is  also  tested  for  the  nonlinear  waves  generated  by  a  simple  source-sink  pair  moving 
below  the  free  surface.  The  aidvantages  of  the  desingularization  is  shown  and  dis¬ 
cussed.  The  desingularized  method  is  then  used  in  the  following  chapters  to  study 
the  wave  problems  mentioned  above. 

Chapter  4,  an  extension  of  the  work  by  Cao,  Schultz  &  Beck  (1991b),  deals  with 
upstream  runaway  solitons  generated  by  a  surface  pressure  or  a  bottom  topography 
moving  near  the  critical  speed  in  two-dimensional  shallow  water.  The  solution  by  the 
desingularized  method  is  compared  to  the  forced  KdV  (fKdV)  solution  (Wu  1987). 
The  difference  between  the  waves  caused  by  the  surface  pressure  and  those  caused  by 
bottom  topography  with  identical  forcing  distribution  in  the  fKdV  model  is  compared 
using  the  fully  nonlinear  model.  The  waves  caused  by  submerged  cylinders,  which 
cannot  be  handled  by  the  fKdV  model,  are  also  studied  using  the  fully  nonlinear 
model. 

Chapter  5,  representing  an  incorporation  of  Cao,  Schultz  &  Beck  (1990),  stud¬ 
ies  the  three-dimensional  waves  generated  by  a  submerged  spheroid  moving  below 
the  free  surface.  The  wave  elevation,  the  pressure  and  forces  on  the  spheroid  are 
calculated  and  compared  to  the  results  by  other  theories  or  methods. 
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In  Chapter  6,  the  mechanism  of  the  generation  of  the  inner-angle  wavepackets  in 
the  wake  of  a  ship  is  studied.  Particulairly,  the  effects  of  the  nonlinear  free  surface 
bound2iry  conditions  and  the  unsteadiness  of  the  waJce  is  examined.  The  wake  is 
studied  with  nonlinear  wave  calculations  using  the  desingularized  method,  linear 
wave  calculations  using  the  time-domain  Green  function  (Liapis  k,  Beck  1986,  King 
1987,  Beck  &  Magee  1990),  and  a  stationary  phase  method.  The  results  of  the 
three  methods  are  compared  and  also  compared  to  the  experimental  data  (Brown  et 
al.  1989).  In  the  nonlinear  calculations,  a  large  computational  domain  on  the  free 
surface  is  required  because  of  our  interest  in  the  far  field,  resulting  in  a  large  number 
of  unknowns  in  the  numerical  computations  (typically  N  >  3000).  Some  numerical 
techniques  for  a  fast  iterative  solution  of  the  discretized  boundary  integral  equation 
are  tested.  These  techniques  include  a  fast  iterative  solver  GMRES  (Saad  1986  aind 
Olson  1989),  a  specially  chosen  precoditioning  matrix,  a  multigrid  method,  as  well 
as  the  vectorization  of  the  computer  program  for  computations  on  supercomputer. 
Finally,  we  make  conclusions  and  recommendations  in  Chapter  7. 


CHAPTER  II 


PROBLEM  FORMULATION  AND  SOLUTION 

PROCEDURE 


In  this  chapter,  the  initial  boundary  value  problem  is  formed  for  nonlinear  prob¬ 
lems  of  waves  generated  by  a  free  surface  pressure  disturbance  or  the  motion  of  a 
rigid  body  moving  imder  the  free  surface.  The  solution  procedure  and  the  numerical 
implementation  are  described.  It  is  assumed  that  the  fluid  is  incompressible  and 
inviscid  and  the  flow  is  irrotational.  This  implies  the  existence  of  a  velocity  potential 
<^(x,  t)  with  the  fluid  velocity  u  given  by  u(x,  t)  =  where  x  =  (x,  y,  z)  is  the 
spatial  location  and  t  is  time.  The  siirface  tension  is  neglected. 

2.1  Initial-Boundary  Value  Problem 

2.1.1  Governing  equation 

The  coordinate  system  used  is  shown  in  Figure  2.1,  in  which  the  z  =  0  plane 
represents  the  mean  position  of  the  free  surface  and  the  z  axis  points  upwards.  The 
problem  under  consideration  is  the  flow  generated  by  the  motion  of  a  rigid  body 
or  singularity  below  the  free  surf^K:e,  or  a  free  surface  pressure  distribution,  or  a 
combination  of  the  two.  The  fluid  domain  D  is  bounded  on  top  by  the  free  surface 
S/,  internally  by  the  hull  of  the  body  5^,  below  by  a  bottom  surface  and  an 
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Figure  2.1:  Problem  definition  and  coordinate  system 
enclosing  contour  S^a  at  infinity. 

Because  of  the  existence  of  the  velocity  potential  the  continuity  equation  re¬ 
quires  that  the  potential  <t>  satisfies  the  Laplace  equation, 

=  0  (in  D).  (2.1) 

It  is  well  known  that  the  momentum  equation  yields  Bernoulli’s  equation  for  the 
irrotational  flow, 

2  =  (inC).  (2.2) 

This  means  that  the  solution  procedures  for  the  velocity  potential  4>  and  the  pressure 
p  are  uncoupled.  Usually,  the  is  first  solved  with  proper  boundary  conditions  and 
initial  conditions.  Then  p  can  be  easily  calculated  using  (2.2). 

2.1.2  Boundary  and  initial  conditions 

Since  the  free  surface  boundary  is  not  known  a  priori,  two  boundary  conditions 
are  required  on  it.  The  kinematic  condition  requires  that  a  fluid  particle  on  the  free 
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surface  will  remain  on  it, 


dr)  d<t>  dr}  d<})d<i>  d<f>_  .  _ 

dt  dx  dx  dy  dy  dz  ^  * 


(2.3) 


where  r}  =  r){x,y,t)  is  the  free  surfaice  elevation.  The  dynamic  condition  requires 
that  the  pressure  on  the  free  surface  equals  the  ambient  pressure  Pa.  neglecting  the 
surface  tension, 

P  =  Paix,  y,  t).  (on  5/)  (2.4) 


The  ambient  pressure  Pa(x,  y,  t)  is  assumed  given  and  can  be  time  and  space  depen¬ 
dent.  Expressing  (2.4)  in  terms  of  (f>  with  the  use  of  (2.2)  on  the  free  surface,  we 
obtain, 

^  (on  5/).  (2.5) 

Eqs.  (2.3)  and  (2.5)  are  the  Eulerian  description  of  the  free  surface  boundary  condi¬ 
tions.  By  re-arranging  Eqs. (2.3)  and  (2.5),  we  obtain  the  Lagrangian  description  of 
the  free  surface  conditions. 


where 


DXj 

Dt 


=  V<f> 


(on  5/) 


D4> 

Dt 


~  -  9^1  + 

P  2 


(on  5/). 


Dt 


V 


(2.6) 

(2.7) 


is  the  substantial  derivative  following  a  fluid  particle  and  Xj  =  (*/,y/,«/)  is  the 
position  vector  of  a  free  surface  particle. 


Depending  on  the  problem  to  be  solved,  one  description  can  be  more  convenient  to 
use  than  the  other.  For  instance,  if  breaking  waves  or  overturning  waves  are  studied, 
the  Lagrangian  description  is  more  convenient  because  the  wave  elevation  becomes 
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multivalued.  The  Eulerian  description  can  be  advantageous  since  the  nodal  spacing 
may  be  kept  uniform  and  free  surface  regriding  is  not  necessary.  In  the  Lagrsmgian 
formulation  the  regriding  may  be  necessary  to  avoid  the  numerical  instabilities  due 
to  Icirgely-uneven  spacing.  We  will  use  the  Lagrangian  description  since  it  has  some 
other  advemtages  in  the  sense  of  the  numerical  wave  computations,  which  will  be 
shown  in  the  following  chapters. 

On  the  body  surface  Sk  and  the  bottom  St,  the  flow  satisfies  the  non-penetration 
condition, 

^  =  (onS,)  (2.8) 

driH 

=  (on5fc),  (2.9) 

OTli 

where  14  and  1^  are  the  velocities  of  body  surface  and  the  bottom  surface,  and  n/, 
and  fit,  are  the  exterior  imit  normal  vectors  to  the  two  surfaces,  respectively. 


At  infinity,  a  proper  far-field  condition  is  required.  For  a  three-dimensional  initial¬ 
boundary  value  problem,  ^  »  0  or  0  implies  that  the  flow  due  to  a  disturbance 

vanishes  at  infinity. 

In  our  time-dependent  problems,  we  study  the  flow  generated  by  disturbances 
starting  from  rest.  Therefore,  for  initial  conditions,  we  require 

^  =  0  (in  i?  for  t  <  0) 

and 


r}  =  T}o{x,y)  {onSj  for  <  <  0), 

where  tJo  is  a  known  initisJ  free  surface  elevation. 


(2.10) 


The  mathematical  formulation  used  here  can  model  the  waves  generated  by  one  of 
four  disturbances  or  any  combination  of  them:  1)  the  motion  or  deformation  of  a  body 
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in  the  water;  2)  the  deformation  of  the  bottom  surface;  3)  the  free  surface  pressure 
distribution;  and  4)  an  initi2il  free  surface  elevation.  For  a  particular  problem,  it  is 
not  necessary  that  all  surfaces  exist.  For  example,  if  the  waves  are  generated  by  the 
free  surface  pressure  disturbance,  surface  5a  does  not  need  to  exist.  Si,  does  not  exist 
if  the  fluid  is  in  deep  water. 

The  mathematical  formulation  is  exact  under  the  potential  flow  assumption.  No 
additional  assumptions  are  made,  such  as  the  linearization  of  the  free  surface  condi¬ 
tions.  The  problem  is  fully  nonlinear  because  of  the  nonlinear  free  surface  conditions. 
The  body  eind  bottom  boundary  conditions  are  exact  although  they  are  linear. 

2.2  Solution  Procedure  for  Nonlinear  Waves 

2.2.1  Solution  procedure 

As  mentioned  in  the  introduction,  two  major  difliculties  encountered  in  the  study 
of  nonlinear  waves  are  that  the  free  surface  boundary  conditions  are  nonlinear  and  the 
free  surface  boundary  is  unknown  a  priori.  There  are  very  few  problems  (usually  with 
simple  boundary  geometry)  for  which  analytic  solutions  can  be  found.  A  numerical 
method  is  necessary  for  practical  computations. 

To  account  for  the  nonlinearity,  a  perturbation  method  can  be  used  that  includes 
higher-order  terms  in  the  series  expansion  of  the  free  surface  boundaury  conditions. 
This  method  has  been  a  powerful  tool  in  the  anal}rtical  study  of  nonlinear  waves 
(Lamb  1929,  Stoker  1957,  Wehausen  &:  Laitone  1960,  amd  Whitham  1974).  The 
perturbation  method  has  also  been  used  in  numerical  calculations  by  researchers,  for 
example,  Nakos  &  Sclavounos  (1990)  in  the  calculation  of  ship  motions  in  waves.  In 
most  cases,  the  boundary  conditions  are  not  satisfied  on  the  read  free  surface  but  on 
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a  flat  mean  free  surface.  An  iterative  procedure  ran  be  used  to  satisfy  the  nonlinear 
free  surface  conditions.  In  the  work  of  Jensen,  Mi  &  Soding  (1986)  and  Jensen, 
Bertram  &  Soding  (1989),  a  linearized  boundary  condition  based  on  the  solution 
from  the  last  iteration  is  satisfied  on  the  deformed  free  surface  determined  from  the 
previous  iteration.  However,  their  methods  are  not  suitable  for  unsteady  waves. 

Longuet- Higgins  &  Cokelet  (1976)  first  introduced  the  mixed  Eulerian-Lagrangian 
method  for  two-dimensional  water  waves  using  fully  nonlinear  free  surface  boundary 
conditions  satisfied  on  the  actual  free  surface  in  numerical  computations.  This  is  a 
time-stepping  procedure  which  requires  two  major  tasks  at  each  time  instant: 

1.  solving  a  boundairy  value  problem,  which  itself  is  linear,  in  an  Eulerian  frame 
to  obtain  the  velocity  on  the  free  surface;  and 

2.  updating  the  positions  of  the  free  surface  particles  and  the  velocity  potential 
on  the  surface  by  integrating  the  nonlinear  Idnematic  and  dynamic  free  surface 
boundary  conditions  with  respect  to  time. 

Variations  of  this  method  have  been  used  for  a  variety  of  two-dimensional  and  three- 
dimensional  nonlinear  free  surface  problems  (Vinje  &  Brevig  1981,  Baker  1982,  Lin 
1984,  Dommermuth  &  Yue  1987,  Kang  1988  and  Cao,  Schultz  &  Beck  1990,  1991a, 
1991b). 

The  boundary  value  problem  at  each  time  step  is. 


t> 

II 

o 

(in  D) 

(2.11) 

O 

II 

■o- 

(on  SfJ 

(2,12) 

9<l>  ^ 

OTlh 

(on  Sh) 

(2.13) 
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-  V  n 
arn 


(on  Si) 


(2.14) 


V^->0 


(as  X  — ♦  oo). 


(2.15) 


where  <i>o  and  S/„  are  known  from  the  previous  time  steps. 


After  solving  the  boundary  value  problem,  the  flow  is  known  for  the  present  time. 
The  velocities  of  the  &ee  surface  particles  can  be  calculated  and  the  free  surface 
position  and  the  potential  <f>o  on  it  can  be  predicted  for  the  next  time  step  by 
integrating  the  free  surface  kinematic  and  dyncimic  conditions  (2.6)  and  (2.7).  The 
pressure  on  the  body  Sh,  and  Sb  can  be  calculated  using  the  Bernoulli’s  equation 


(2.2), 


where 


=  -  (j  V  *  -  V)  ■ 


(on  Sk  or  Si), 


(2.16) 


is  the  substantial  derivative  of  the  potential  at  a  fixed  point  on  Sh  or  S\,  following 
the  motion  of  that  point  and  V  is  the  velocity  of  that  point  (i.e.  V  =  Vh  or  H).  The 
second  form  of  (2.16)  is  more  appropriate  to  use  because  the  nodal  points  must  move 
to  follow  the  boundary.  The  forces  and  the  moments  on  the  body  can  be  calculated 
by  integrating  the  pressure  over  the  the  body  surface, 


(2.17) 


M  =  J  J  — p(r  X  nk)ds. 


(2.18) 
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where  r  is  the  position  vector  of  a  body  surface  point  to  a  reference  point  (usually 
the  centroid  of  the  body). 

2.2.2  Numerical  implementation 

The  boundary  value  problem  (2.11  —  2.15)  has  to  be  solved  by  a  numerical 
method.  Many  methods,  such  as  finite  difference  methods,  finite  element  methods 
and  boundary  integral  methods,  are  available.  A  review  of  numerical  methods  for 
free  surface  flows  is  given  in  Yeung  (1982).  It  is  very  important  to  choose  an  ef¬ 
fective  solution  method  for  the  boimdary  value  problem  since  its  solution  requires 
most  of  the  computation  time.  A  boundary  integral  method  is  powerful  because  it 
reformulates  the  origin2d  problem  as  an  integral  equation  on  the  domain  boimdary 
thus  reducing  the  computational  domain  by  one  dimension.  We  use  a  desingularized 
boundary  integral  method  to  solve  (2.11  —  2.15).  Although  the  method  is  not  as 
widely  used  as  conventional  singular  boundary  integral  methods  because  of  lack  of 
a  uniqueness  and  completeness  theorem,  it  has  many  potential  advantages  over  the 
singular  methods  and  has  received  more  attention  recently.  In  Chapter  3,  the  desin¬ 
gularized  method  is  re-examined  and  tested  with  some  simple  problems.  It  is  shown 
that  the  desingularized  method  is  faster,  easier  to  implement  (especially  easy  to  take 
advantage  of  vector  and  parallel  computers)  and  can  be  more  accurate  for  many 
problems.  Therefore,  this  method  is  used  in  the  nonlinear  wave  problems  studied  in 
Chapters  4,  5  and  6. 

The  time  integration  of  the  free  surface  boundary  conditions  (2.6)  and  (2.7), 
which  are  first-order  ordinary  differential  equations  with  respect  to  time,  are  inte¬ 
grated  with  a  widely  used  fourth-order  Runge-Kutta-Fehlberg  method.  This  method 
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has  been  shown  to  be  very  accurate  and  stable.  The  subroutine  can  be  found  in  ODE- 
PACK.  An  initial  time  increment  is  set,  but  is  modified  by  the  Runge-Kutta-Fehlberg 
subroutine  where  appropriate  to  satisfy  an  error  tolerance. 


CHAPTER  III 


DESINGULARIZED  BOUNDARY  INTEGRAL 

METHOD 


A  boundary  integral  equation  method  employs  a  fundamental  solution,  which 
satisfies  the  differential  equation  (and  possibly  part  of  the  boundary  conditions), 
to  reformulate  the  original  boundary  value  problem  as  an  integral  equation  on  the 
boundary.  In  conventional  boundary  integral  formulations,  singularities  of  the  funda¬ 
mental  solution  are  placed  on  the  domain  boundary.  This  requires  special  evaluation 
of  singular  integrands,  which  can  result  in  costly  numerical  calculations.  In  the 
time-dependent  procedure  for  nonlinear  free  surface  problems,  the  boundary  integral 
equation  needs  to  be  solved  at  each  time  step.  Since  most  of  the  computation  time 
is  devoted  to  the  boundary  integral  problem,  an  effective  solution  method  is  critical 
in  the  time-stepping  procedure. 

When  the  singularity  of  the  fundamental  solution  is  placed  away  from  the  boimd- 
ary  and  outside  the  domain  of  the  problem,  a  desingularized  boundary  integral  equa¬ 
tion  is  obtained.  Similar  to  the  singular  formulation  of  boundary  integral  equations, 
the  non-singular  formulation  also  has  two  versions:  the  direct  method  and  the  in¬ 
direct  method.  In  the  direct  method.  Green’s  second  identity  is  used  to  derive  the 
boundary  integral  equation,  and  the  solution  of  the  problem  is  obtained  directly 
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by  solving  the  boundary  integral  equation.  In  the  indirect  method,  the  solution  is 
constructed  by  a  singularity  distribution  outside  the  problem  domain  and  a  bound¬ 
ary  integral  equation  for  the  strength  of  the  singularities  is  formulated  to  satisfy 
the  boundary  condition.  After  the  strength  is  determined,  the  flow  can  be  easily 
calculated. 

One  immediate  advantage  of  the  desingularization  is  that  simple  quadrature  can 
be  used  to  evaluate  the  integrals  since  the  integrals  are  non-singular,  thus  reducing 
the  computational  time  in  obtaining  the  algebraic  system  representing  the  discretized 
boundairy  integral  equation.  It  has  also  been  shown  that  the  desingularized  boundary 
integral  equation  can  give  more  accurate  solutions  for  a  given  tnmcation  in  some 
problems  (Webster  1975,  Schultz  &  Hong  1989). 

The  first  use  of  a  desingularization  method  is  the  classical  work  by  von  Karman 
(1930)  that  determines  the  flow  about  an  axisymmetric  body  using  an  axial  source 
distribution.  The  strength  of  the  source  distribution  is  determined  by  the  kinematic 
boundary  condition  on  the  body  surface.  Kupradze  (1967)  proposes  a  non-singular 
formulation  of  the  direct  boimdary  integral  equation  for  the  exterior  Dirichlet  prob¬ 
lem  with  an  auxiliary  control  surface  outside  the  problem  domain.  Heise  (1978) 
studies  some  numerical  properties  of  integral  equations  in  which  the  singular  points 
are  on  an  auxiliary  boundary  outside  the  solution  domain  for  plane  elastostatic  prob¬ 
lems.  By  applying  Green’s  theorem  to  the  solution  and  a  simple  linear  wave  source 
with  the  singular  point  lying  inside  the  body  (i.e.,  outside  the  problem  domain)  and 
using  a  bilinear  expansion  of  the  source,  Martin  (1981)  obtains  the  so-called  “null- 
field  equations  for  the  water  wave  radiation  problenos”.  Johnston  k  Fairweather 
(1984)  and  Han  k  Olson  (1987)  use  an  adaptive  method  in  which  the  singularities 
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axe  located  outside  the  domain  and  allowed  to  move  as  part  of  the  solution  process. 
This  adaptive  method  requires  considerably  fewer  singularities  than  the  number  of 
boundary  nodes,  but  it  results  in  a  system  of  nonlinear  algebraic  equations  for  both 
the  strength  and  the  location  of  the  singularities. 

Webster  (1975)  gives  a  discussion  on  the  numerical  properties  of  the  desingular- 
ization  technique  for  the  external  potential  flow  of  an  arbitrary,  three-dimensional 
smooth  body.  Triangular  panels  of  a  simple  source  distribution  are  placed  some¬ 
what  inside  the  surface  of  body.  The  integration  is  performed  analytically  for  each 
panel  using  a  bilinear  distribution  of  the  singularity  strength.  These  integrations 
require  evaluation  of  (logarithmic  and  arctangent)  transcendental  functions.  From 
the  numerical  results,  Webster  concludes  that  “submergence  of  the  singularity  sheet 
below  the  surface  of  the  body  appears  to  improve  greatly  the  accuracy,  as  long  as  the 
sheet  is  not  submerged  too  far”.  This  indicates  that  for  a  certain  discretization  of 
the  body  surface,  a  more  accurate  solution  may  be  obtained  using  the  non-singular 
formulation  than  the  singular  formulation.  Although  Webster  did  not  attempt  it, 
a  simple  numericaJ  quadrature  can  be  used  if  the  desingularization  distance  is  suf¬ 
ficiently  large.  This  will  significantly  reduce  the  computational  effort  to  obtain  the 
algebraic  system. 

Although  the  non-singular  formulation  of  the  indirect  method  has  grown  popular 
recently,  few  studies  on  the  non-singular  direct  method  have  been  reported,  espe¬ 
cially  for  three-dimensional  problems.  Schultz  &  Hong  (1989)  obtain  a  non-singular 
direct  formulation  by  moving  the  singularity  in  the  Cauchy’s  integral  away  from  the 
boundary  in  two-dimensional  potential  problems.  Their  results  show  the  advantages 
of  the  non-singular  formulation.  They  also  use  an  overdetermined  system  combin- 
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ing  the  real  and  imaginary  parts  of  Cauchy’s  theorem  or  using  “extra”  evaluation 
points  away  from  the  boundary  contour.  This  overdetermined  system  can  exhibit 
higher-order  convergence  than  the  determined  system  from  either  the  real  part  or 
the  imaginary  part  of  Cauchy’s  theorem.  Few  investigations  on  applications  with 
boundaries  extending  to  infinity,  espedadly  for  water  wave  problems,  have  been  re¬ 
ported.  The  work  of  Jensen,  Mi  fz  Soding  (1986)  solving  the  steady  nonlinear  ship 
wave  problem  using  a  Rankine  source  distribution  above  the  free  surface  can  be 
considered  as  an  application  of  the  indirect  desingularized  method. 

The  direct  and  indirect  desingularized  methods  are  formulated  in  the  following 
sections.  The  effect  of  desingularization  distance  is  examined  for  a  simple  problem 
with  an  infinitive  and  flat  boundary.  The  convergence  of  the  method  with  respect 
to  mesh  size  and  computational  time  are  compared  to  find  the  “optimum”  desingu¬ 
larization  distance  (will  be  discussed  in  more  detsul  later).  This  “optimum”  distance 
will  then  be  used  for  the  simulations  of  unsteady  nonlinear  waves  in  the  following 
chapters. 

3.1  Desingularized  Integral  Equations 

The  Laplace  equation  is  the  governing  equation  in  a  domain  D  bounded  by  F.  A 
Dirichlet  condition  is  imposed  on  part  of  the  boundaury  F ^  amd  a  Neumann  condition 
on  the  remaining  boundary  F„.  Therefore,  one  has  the  following  mixed  boundary 
value  problem: 


> 

li 

o 

(in  D), 

(3.1) 

II 

o 

(on  Fj), 

(3.2) 

II 

X 

(on  r„). 

(3.3) 

where  and  x  are  known  functions.  If  the  boundary  extends  to  infinity,  a  far-held 
condition  is  required.  The  desingularized  boundary  integral  method  separates  the 
integration  and  control  (i.e.,  evaluation)  surfaces,  one  of  which  is  the  boundary  of 
the  problem.  In  the  direct  method,  the  boundary  of  the  problem  is  the  integration 
surface,  while  in  the  indirect  method,  the  boundary  is  the  control  surface. 


3.1.1  Direct  method 


The  boundary  integral  equation  in  the  direct  method  is  derived  &om  Green’s 
second  identity: 

III  -  M*)  dD-ll  =  0,  (3.4) 

which  holds  for  any  two  functions  with  second  derivatives  continuous  in  D  and  the 
boundary  F.  If  4  is  the  solution  and  tl>  is  chosen  as 


(3.5) 


with  its  evaluation  point  Xj,  outside  D  and  integration  point  X,  on  F,  the  volume 
integral  in  (3.4)  vanishes.  Equation  (3.4)  then  becomes 


Applying  the  boundary  conditions  (3.2)  and  (3.3)  to  (3.6)  and  moving  the  known 
quantities  to  the  right-hand  side  gives: 


d<t>{x,) 


Xj,-X, 


dn 


dr 
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The  kernels  are  non-singular  since  Xp  and  X,  never  coincide.  The  int^ral  equation 
(3.7)  is  solved  for  <f>  on  Fn  and  d4>ldn  on  F^. 

The  fundamental  solution  |  can  easily  be  replaced  by  other  higher-order 

singularities  (dipole,  etc.)  with  little  additional  computational  effort  since  the  inte¬ 
grals  are  non- singular  and  can  be  evaluated  numerically. 
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3.2  Uniqueness  of  Desingularized  Integral  Equation  and 
Convergence  of  the  Numerical  Methods 

The  existence  of  the  solution  to  the  direct  version  of  the  desingularized  boundary 
integral  equation  is  obvious  as  long  as  the  original  boundary  value  problem  is  well- 
posed  and  has  a  unique  solution  according  to  the  Green’s  theorem. 

Kupradze  (1967)  gives  a  proof  of  the  uniqueness  of  the  solution  to  the  direct  ver¬ 
sion  of  the  desingularized  boundary  integral  equation  for  exterior  Dirichlet  problems 
£Lnd  proposes  an  approximate  solution  method  in  which  the  solution  is  expressed  as 

(3.11) 

where 

'^kiXp)  —  T-^  ^  r  (^  =  1>2,...)  (3.12) 

is  a  linearly  independent  and  complete  set.  Xp  is  the  point  on  the  real  boundary  and 
(k  =  1, 2, 3, ...)  are  points  covering  the  auxiliary  boundary  everywhere  “densely” . 
Substituting  (3.11)  into  the  desigularized  equation  gives  a  system  of  linear  algebraic 
equations  for  Ck  {k  =  1,7V).  This,  in  fact,  is  a  kind  of  collocation  method  where 
is  the  collocation  point.  Kupradze  proves  the  convergence  of  the  method.  However, 
this  method  may  not  be  very  effective  in  practice  since  each  element  of  the  coefficient 
matrix  and  the  right  hand  side  of  the  linear  system  requires  an  integration  over  the 
whole  real  boundary  that  will  still  be  very  time  consuming  although  the  integrands 
are  non-singular.  Also,  in  his  method,  the  auxiliary  boundary  does  not  change  as  N 
changes,  which  may  cause  numerical  instabilities  as  N  becomes  large. 

A  simple  collocation  method  is  used  here  to  solve  the  integral  equations.  In 
the  direct  method  the  Green’s  theorem  is  evaluated  at  chosen  collocation  points  on 
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the  auxiliary  boundary,  while  in  the  indirect  method  the  boundary  conditions  are 
satisfied  at  the  collocation  points  on  the  real  boundary.  The  integration  surface  (the 
real  boundary  in  the  direct  method  or  the  auxiliary  boundsury  in  the  indirect  method) 
is  divided  into  panels,  a.nd  the  integration  is  carried  out  for  each  panel  assuming  some 
kind  of  distribution  of  the  solution  or  source  strength.  Therefore,  each  element  in 
the  coefficient  matrix  of  the  resulting  linear  syston  requires  integrations  over  a  small 
number  of  the  nearby  panels.  This  will  reduce  the  computation  time  signific2Lntly. 

The  convergence  of  a  numerical  method  requires  that  the  numerical  solution 
approach  the  exact  solution  as  the  computational  grid  is  refined.  To  ensure  the  con¬ 
vergence,  the  non-dimensional  desingularization  distance  of  a  singulsir  point  (outside 
the  domain  and  normal  to  the  boundary)  is  proposed  to  be  given  by 

L4  =  (3.13) 

where  I4  is  a  parameter  independent  of  the  discretization  (/^  will  be  termed  the 
desingularization  factor)  which  reflects  how  far  the  integral  equation  is  desigularized. 
Dm  is  the  local  mesh  size  (usually  the  square  root  of  the  local  mesh  area  in  a  three- 
dimensional  problem)  and  1/  is  a  parameter  that  affects  the  accuracy  of  the  munerical 
integration. 

The  use  of  (3.13)  implies  that  the  desingularization  distance  is  related  to  the 
local  mesh  size.  As  the  mesh  becomes  finer,  the  singular  point  becomes  closer  to 
the  boundary  for  1/  >  0  .  In  the  limit,  the  non-singular  formulation  is  consistent 
with  the  singxilar  formulation  although  the  kernels  do  not  become  singular  and  the 
desigularized  integr2d  equations  are  expected  to  have  simileu*  properties  ais  the  singular 
boundary  integral  equations.  Numerically,  the  use  of  the  desingularization  distance 
(3.13)  with  1/  >  0  is  also  advantageous  over  the  case  when  u  =  0  since  the  condition 
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number  of  the  coefficient  matrix  of  the  desingularized  equations  will  deteriorate  less 
rapidly  as  the  mesh  is  refined. 

Desingularization  increases  the  condition  number  of  the  coefficient  matrix  for  a 
given  truncation.  Webster  (1975)  shows  that  non-singular  formiilations  are  not  sig¬ 
nificantly  less  well  conditioned  than  singular  formulations  if  the  desingularization 
is  not  too  large.  In  most  three-dimensional  computations,  the  large  number  of  un¬ 
knowns  requires  an  iterative  solution  of  the  linear  system.  As  the  desingularization 
increases,  the  condition  number  also  in'^reases.  Thus,  an  increase  in  the  number  of 
iterations  in  solving  the  algebraic  system  is  expected.  There  appears  to  be  an  “op¬ 
timum”  desingularization  distance:  if  the  singularity  is  too  f^  from  the  boundary, 
the  linear  system  will  be  poorly  conditioned  and  numerical  instabilities  occur;  if  the 
singularity  is  too  close  to  the  boundary,  numerical  integration  of  the  increasingly 
singular  kernel  is  suspect  and  the  solution  may  not  be  accurate.  The  “optimum” 
and  proper  u  will  be  determined  by  numerical  tests  for  simple  problems. 

3.3  Simple  Test  Examples 

To  test  the  convergence  of  the  numerical  integration  with  respect  to  the  mesh 
refinement,  the  potential  due  to  a  known  constant  normal  dipole  distribution  within 
a  square,  flat  surface  (—1  <  i  <  1,  —1  <  y  <  1,  r  =  0)  is  considered.  The  potential  is 
evaluated  at  a  point  above  the  center  of  the  surface  at  a  distance  given  by  (3.13)  with 
Id  =  1.  The  surface  is  divided  into  square  panels  and  a  2  x  2  Gaussian  quadrature  is 
used  for  each  panel.  A  third-order  convergence  should  be  found  assuming  the  inte¬ 
grand  is  non-singular  and  is  independent  of  the  mesh  (i.e.  i/  =  0).  For  i/  >  0,  (3.13) 
makes  the  integrand  depend  on  the  mesh  size,  and  hence  the  third-order  convergence 
is  not  assured.  Fig.  3.1  shows  the  convergence  of  the  numerical  integration  for.  three 
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N 

Figiire  3.1:  Convergence  of  numerical  integration 


values  of  vx  0,  1/2  and  1.  As  can  be  seen,  the  numerical  integration  error  increases 
with  V.  However,  u  should  be  greater  than  0  for  the  convergence  of  the  solution  of 
the  integral  equation.  Therefore,  an  appropriate  u  lies  in  between  0  and  1.  i/  =  1/2 
is  chosen  for  the  subsequent  calculations. 

The  numerical  performance  of  the  desingulairized  boundary  integral  methods  is 
tested  for  a  problem  in  which  the  potential  is  generated  by  a  dipole  below  an  infinitive 
flat  plane  with  ^  =  0.  This  simple  problem  has  an  exaict  solution  formed  by  the  dipole 
aind  its  negative  image  about  the  flat  plane.  This  problem  represents  the  solution  to 
the  first  time  step  of  nonlineair  waves  caused  by  am  impulsive  disturbance  of  a  dipole 
under  water.  The  method  is  then  tested  for  simple  time-dependent  nonlinear  waves 
caused  by  a  moving  source-sink  pair  under  the  free  surface. 
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The  direct  and  indirect  methods  axe  compared  in  this  example.  In  the  direct 
method,  the  boundary  is  divided  into  panels,  and  a  solution  is  sought  at  the  nodes. 
Since  the  kernels  are  non-singular,  the  panel  integrals  can  be  evaluated  using  Gaus¬ 
sian  quadrature,  where  the  solution  is  assumed  to  be  bilinear  over  each  panel.  In  the 
indirect  method,  the  integrals  of  the  singularity  distribution  (3.8)  is  replaced  by  a 
summation  of  concentrated  singularities  without  an  apparent  loss  of  accuracy  when 
the  desingularization  distance  is  appropriate.  Therefore,  the  computations  will  be 
less  complex  and  time  consuming  than  in  the  direct  method. 

In  the  results  that  follow,  the  Generalized  Minimal  Residual  Method  (GMRES), 
Saad  &  Schultz  (1986)  and  Olson  (1989),  is  used  to  iteratively  solve  the  system  of 
equations.  It  is  found  that  this  method  to  be  generally  twice  as  fast  as  a  standard 
conjugate  gradient  routine  for  nons}rmmetric  matrices.  A  sufficiently  small  value  of 
the  convergence  tolerance  is  chosen  so  that  it  introduces  negligible  error  in  compari¬ 
son  to  that  introduced  by  truncation.  The  computations  are  carried  out  using  double 
precision  on  an  Apollo  DNIOOOO  workstation. 

3.3.1  A  dipole  below  a  ^  =  0  infinitive  fiat  plane 

In  this  example,  a  Dirichlet  condition,  ^  =  0,  is  imposed  on  the  z  =  0  plane. 
A  dipole  of  unit  strength  is  located  at  .^o(0,0,  — 1),  and  the  direction  of  the  dipole 
coincides  with  the  x  axis.  The  normsd  derivative  d(^fdn  is  sought  on  x  =  0.  The 
problem  is  symmetric  about  the  y  =  0  plane.  The  computations  take  advantage  of 
this  symmetry. 

Since  the  potential  at  the  dipole  location  is  singular  for  this  problem,  the  solution 
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is  decomposed  into  a  regular  part  and  a  singular  part: 


(3.14) 


For  the  direct  method,  we  apply  (3.7)  for  4  recognize  the  r„  is  absent  and 
far-held  closure  does  not  contribute.  Substituting  back  for  <t>  gives 

d<t> 
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(3.15) 


/r-  dn  |je,  -  ;e,|  -  X 

where  Xp  =  (xp,yp,Zp)  is  the  control  point  outside  D,  and  Xg  the  integration  point 
on  Fj  (z  =  0).  The  computational  domain,  the  truncation  of  F^  by  x  =  ±Xo  and 
y  =  ±i2oo,  is  discretized  into  a  uniform  square  mesh  defined  by  N  node  points. 
The  mesh  configuration  is  shown  in  Fig.  3.2.  The  integrations  assume  a  bilinear 
distribution  of  d(f>fdn  and  use  a  Gaussian  quadrature.  The  control  points  Xp  are 
placed  directly  above  the  node  points  at  a  distance  =  /j(AA^)^/^,  where  Ay  is 
the  side  of  a  square  .panel. 

For  the  indirect  method,  we  construct  the  potential  using  (3.9),  with  the  integra¬ 
tion  replaced  by  summation  as  discussed  previously,  supplemented  with  the  singulzu’ 
dipole  term 


<t>{X)  = 


-1 
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|3  +  S  I  ^ 


<y. 


(3.16) 


4.|;e-xr  ’ 

where  Xn  are  the  source  points  above  z  =  0.  The  boundary  condition  ^  =  0  on 
z  =  0  results  in 

N 


^  <Ti  1 


where  X  =  {xc,yeiZc)  is  the  control  point  on  z  =  0. 


(3.17) 


The  exact  solution  for  <f>{X)  is 
^(x,y,z)=  ^  ^ 


Air 


+  1,2  +  (r  +  1)2) (l2  +  5,2  +  _  1)2) 
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(3.18) 
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Figure  3.2:  Schematic  diagram  of  a  dipole  below  a  ^  =  0  infinitive  flat  plane 


The  exact  normal  derivative  of  <f>  on  the  z  =  0  plane  then  is 

_ 3 _ X _ 


(3.19) 


The  numerical  solution  for  the  normal  derivative  of  is  comp2ired  to  (3.19)  by  the 


RMS  error  measurement  defined  by 


(3.20) 


Fig.  3.3  shows  the  RMS  error  as  a  function  of  for  three  truncations.  In  the 
direct  method,  in  which  a  2  x  2  Gauss  quadrature  is  used  to  evaluate  the  influence 
function,  the  RMS  error  increases  as  increases.  In  the  indirect  method,  the  solution 
breaks  down  Ij  =  0  because  the  integration  and  control  points  coincide.  However, 
as  Ij  increases,  the  RMS  error  decreases  rapidly  and  remains  rather  insensitive  to 
Id  for  quite  a  large  range.  Eventually,  Ij  becomes  too  large  for  the  given  truncation 
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to  represent  the  second  term  in  (3.18)  well,  and  the  error  starts  to  increase.  It  is 
interesting  that  desingularization  in  the  indirect  method  greatly  reduces  error  and 
gives  smaller  error  than  in  the  direct  method  for  a  large  range  of  One  of  the 
reasons  is  that  solution  of  many  problems  can  be  represented  accurately  by  a  finite 
number  of  some  singularities  located  outside  the  problem  domain.  In  this  example, 
one  negative  image  dipole  above  the  a  =  0  plane  results  in  the  exact  solution.  The 
direct  method  seems  to  be  less  sensitive  to  the  desingularization. 

The  effect  of  different  quadratures  to  evaluate  the  coeflBcient  matrix  for  the  direct 
method  is  given  in  Fig.  3.4.  As  expected,  for  small  /j,  the  analjrtic  integration  by 
the  approach  given  in  Newman  (1986)  produce  smaller  RMS  error  than  the  Gaussian 
quadrature  because  the  integration  error  dominates  the  accuracy  of  the  solution.  As 
increases,  the  Gaussian  quadrature  integrates  the  smoother  integrands  accurately, 
and  the  results  merge  to  those  using  the  analytic  integration.  It  is  fortuitous  that 
in  the  intermediate  range  of  U,  the  results  using  the  Gaussian  quadrature  show  less 
error  than  those  using  the  2uialytic  integration — the  discretization  and  quadrature 
errors  seem  to  pairtially  cancel  each  other.  A  more  accurate  3x3  Gaussiem  quadrature 
brings  the  results  closer  to  those  using  the  analytic  integration,  as  expected.  For  the 
singular  direct  method  (/,j  =  0  in  Fig.  3.4),  the  analytic  integration  gives  the  least 
error  as  expected. 

The  effect  of  desingularization  on  the  condition  number  of  the  algebraic  system  is 
shown  in  Fig.  3.5.  As  expected,  the  condition  number  increases  with  Id-  However,  a 
poorly-conditioned  system  does  not  necessarily  imply  an  inaccurate  solution  as  long 
as  Id  is  not  too  large.  In  our  example,  minimal  error  occurs  around  Id  =  3.0  in  the 
indirect  method,  where  the  condition  number  is  quite  large. 
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Figure  3.3:  Effect  of  desingularization  (Hoo  =  6.6667): 
o:N  =  231;  a:N  =  496;  f:N  =  861; 

- direct  method; - indirect  method 


0  0.5  1  1.5  2 


Figure  3.4:  Comparisons  between  Gaussian  quadrature  and  anal)d;ic  integration 
(Roo  =  6.6667,  N  =  231) 
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Figure  3.5:  Effect  of  desingularization  on  condition  number 
(Rco  =  6.6667,  N  =  231) 

Fig.  3.6  compares  the  ratios  of  the  CPU  time  by  the  GMRES  to  that  by  an  LU 
algorithm  for  the  direct  method  with  231  unknowns  as  a  function  of  The  LU 
algorithm  takes  about  3.5  seconds  to  solve  the  matrix  equation.  An  error  tolerance  c 
for  the  least  squares  residual  of  the  equations  needs  to  be  specified  when  the  GMRES 
IS  used.  We  find  that  c  =  10“®  is  sufiSciently  small  in  all  the  computations  performed 
in  this  example.  With  the  use  of  this  tolerance,  the  CPU  time  by  the  GMRES  is 
less  than  that  by  the  LU  algorithm  for  U  <  1.5.  The  analytic  integration  requires 
about  22  seconds  to  set  up  the  matrix  while  the  2  x  2  Gaussian  quadrature  requires 
only  about  6.4  seconds,  for  a  savings  of  about  70  percent.  The  CPU  time  for  solving 
the  matrix  equation  by  the  GMRES  varies  from  3.0  seconds  to  4.0  seconds  as  U 
changes  from  0.8  to  3.0.  Even  in  the  worst  case,  U  =  3.0,  we  have  a  total  reduction 
m  the  CPU  time  of  about  60  percent  over  the  singular  formulation  (with  /j  =  0  and 
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Figure  3.6:  Effect  of  desingularization  on  computation  time 
(Eoo  =  6.6667,  N  =  231) 


analytic  integration). 

For  the  indirect  method  with  the  same  truncation,  the  CPU  time  for  the  matrix 
set  up  is  only  about  1.3  seconds,  and  the  CPU  time  for  the  solution  of  the  matrix 
equation  by  the  GMRES  varies  from  1.7  to  4.0  seconds.  The  total  CPU  time  is  re¬ 
duced  by  about  80  percent  using  the  non-singular  indirect  method  for  this  truncation 
compared  with  the  singular  direct  method.  When  Id  is  small  (e.g.,  <  1.5),  a  larger  c 
does  not  result  in  a  significant  difference  in  the  RMS  error  (Fig.  3.7),  but  the  CPU 
time  can  be  significantly  reduced  (Fig.  3.8). 

Fig.  3.9  shows  the  effect  of  the  tnmcation  of  the  infinitive  plane.  The  computa¬ 
tional  domain  is  extended  by  adding  uniformly  sized  meshes.  As  expected,  both  the 
direct  and  indirect  methods  converge  quadratically  with  respect  to  the  extent  of  the 
computational  domain  for.  this  problem.  Both  methods  also  converge  algebraically 


Figure  3.7:  Effect  of  iteration  tolerance  on  RMS  error  {Roo  =  6.6667,  N  =  S 
/<i  =  1,  indirect  method) 


Figure  3.8:  Effect  of  iteration  tolerance  on  computation  time  {R^  =  6.6667,  N 
231,  /rf  =  1,  indirect  method) 
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(approximately  linearly)  with  respect  to  iV,  Fig.  3.10,  in  which  the  computational 
domain  remains  unchanged  while  the  mesh  size  varies. 

We  also  use  a  surface  distribution  of  sources  above  the  2  =  0  plane  in  the  indi¬ 
rect  method  for  this  problem.  This  is  similar  to  Webster’s  method  (Webster  1975) 
except  that  the  integr2ds  over  the  source  surface  are  carried  out  by  a  2  x  2  Gaussian 
quadrature  in  our  case.  As  seen  in  Fig.  3.11,  the  results  axe  better  than  those  using 
concentrated  sources  for  a  large  range  of  (from  0.1  to  3.0).  This  cam  be  expected 
since  the  use  of  concentrated  sources  implies  using  a  rougher  quadrature  for  the  in- 
tegrads  over  the  source  surface.  However,  we  found  that  the  condition  number  of  the 
resulting  system  of  equations  using  a  surface  distribution  of  sources  is  greater  than 
that  using  concentrated  sources  by  0(10^)  (not  shown  here).  Therefore,  the  system 
is  no  longer  well  conditioned  for  /</  greater  than  3.0  while  the  condition  of  the  system 
using  concentrated  sources  is  still  acceptable  with  a  double  precision  computation. 
In  addition,  the  distributed  source  requires  about  as  4  times  the  CPU  time  as  the 
concentrated  sources  in  the  matrix  setup. 

Fig.  3.12  shows  the  error  distribution  along  the  x  axis  for  both  direct  and  indirect 
methods.  Larger  oscillations  in  the  solution  by  the  direct  method  are  observed  at 
the  edge  of  the  computational  domain.  This  may  be  due  to  the  neglecting  of  the 
contribution  from  the  integrals  over  the  far-field  closure  in  the  direct  method,  which 
is  very  likely  to  result  in  a  strong  global  effect  on  the  accuracy  of  the  boundary 
integral  equation  itself.  In  the  indirect  method  no  integr£ds  are  required  over  the  real 
boundary;  the  boundary  condition  is  enforced  within  the  computational  dom£un,  and 
the  singularities  individually  satisfy  the  far-field  condition.  Thus,  one  may  expect 
that  the  effect  due  to  the  failure  in  satisfying  the  boundary  condition  outside  the 
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Figure  3.11:  Compaxison  of  surface  distribution  and  isolated  sources  for  indirect 
method  (R^  =  6.667,  N  =  231) 

computational  domain  to  be  smaller  than  the  effect  due  to  the  neglecting  of  the 
far-field  closure.  Overall,  we  find  the  indirect  method  to  be  better  than  the  direct 
method  for  the  problems  with  infinitive  boundaries. 

3.3.2  Waves  due  to  a  source-sink  pair  moving  below  the  free  surface 

In  this  example,  the  waves  are  generated  by  the  motion  of  a  submerged  source- 
sink  pair  below  the  free  surface.  The  problem  is  nondimensionalized  by  making 
the  submergence  depth  of  the  source-sink  pair,  the  gravity  acceleration  and  the 
fluid  density  unity.  The  source-sink  pair  moves  in  the  x  direction  with  the  velocity 
^  ~  -  e  starting  from  rest  and  gradually  increasing  to  a  steady  value  F^. 

The  strength  of  the  pair  is  given  by  <T{t)  =  <t,(1  -  c’'**)  which  also  increases  from 
zero  to  the  final  value  (To  gradually.  The  distance  between  the  source  and  sink  is  set 
equal  to  0.1.  The  Froude  number  is  unity.  The  midpoint  between  the  pair  is  initiaUy 
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Figure  3.12:  Error  distribution  along  z-axis  (iJoo  =  6.6667,  N  =  231,  Id  =  1) 
located  at  point  (5,0,— 1). 

We  use  the  indirect  method  with  concentrated  sources  because  of  its  computa¬ 
tional  advantages  and  simplicity.  An  additional  advantage  of  the  indirect  method 
is  that  the  velocities  can  be  calculated  directly  once  the  strength  of  the  singularity 
distribution  is  known,  while  the  direct  method  requires  numerical  differentiation  to 
obtain  the  tangential  components  of  the  velocity.  The  potenticil  is  expressed  as  a 
sum  of  1)  the  potential  due  to  the  submerged  source-sink  pair  in  a  infinitive  fluid,  2) 
the  image  of  the  submerged  pair  above  the  undisturbed  free  surface,  and  3)  a  sum 
of  N  concentrated  sources  of  unknown  strength  above  the  free  surface.  The  fax-field 
condition  is  better  satisfied  when  the  image  of  the  pair  is  used  in  the  construction  of 
the  solution.  The  discretized  integral  equation  for  the  unknowns  <Ti  is 
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where  4>o{.Xf)  is  known  from  previous  time.  J^touree  and  Xsink  are  the  locations  of  the 
source  and  sink  and  Xl^^^  and  are  the  locations  of  their  images.  After  the  <Ti 
are  solved,  the  velocities  of  the  fluid  particles  on  the  free  surface  can  be  calculated 
A  fourth-order  Runge-Kutta-Fehlberg  method  is  used  in  the  nonlinear  free  surface 
integration. 


In  the  results  shown  below,  the  initial  free  surface  mesh  grid  extends  from  0  < 
I  <  20  and  0  <  y  <  7.5  and  has  41  x  16  node  points.  The  grid  spacing  is  uniform  in 
the  X  direction  but  increases  by  10%  in  the  y  direction.  The  node  points  serve  as  the 
material  points.  The  concentrated  sources  are  placed  approximately  perpendicular 
from  the  node  points  at  a  distance  determined  by  (3.13)  with  /j  =  1.0  and  u  =  0.5. 

The  method  is  first  applied  to  waves  generated  by  a  sufficiently  small  disturbance 
such  that  hnear  wave  theory  is  a  good  approximation.  The  results  of  the  present 
method  using  fully  nonlinear  free  surface  conditions  are  compared  to  an  “exact”  solu¬ 
tion  computed  from  a  time-dependent  Green  function  for  a  Kelvin  wave  source  that 
satisfies  the  linearized  free  surface  condition  (King,  Beck  &  Magee,  1988).  Fig.  3.13 
shows  the  comparison  of  the  nonlinear  wave  elevation  along  the  symmetry  plane  to 
the  linear  wave  elevation  at  t  =  10  for  a  weak  disturbance  (o-^  =  0.05).  The  nonlinear 
and  linear  results  agree  very  well.  Independent  computations  using;  a)  a  smaller  com¬ 
putational  domain  (with  the  same  mesh  spacing  within  0  <  i  <  15  and  0  <  y  <  7.5), 
b)  finer  mesh  grids  (81  x  16  and  41  x  31  with  the  same  computational  domain  size), 
and  c)  doubling  the  time  increment,  result  in  negligible  difference  for  the  nonlinear 
calculation.  This  indicates  that  even  for  the  small  disturbance  example  studied  here, 
the  differences  in  Fig.  3.13  are  primarily  due  to  nonlinear  effects.  Fig.  3.14  shows 
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the  results  for  a  stronger  disturbance  ((To  =  0.75),  showing  the  significant  nonlinear 
effects  of  the  free  surface  conditions,  especially  at  the  trough. 

To  study  the  effects  (maiinly  the  non-physi<»l  reflection)  due  to  the  truncation 
of  the  free  surface,  we  compare  the  results  using  two  differently  sized  computational 
domains  mentioned  above  with  identical  grid  spacing  for  the  weak  disturbance  case. 
Fig.  3.15  shows  the  results  using  the  two  domains.  Little  difference  is  seen  in  the 
results  of  the  two  computations  until  the  front  wave  hits  the  boundary  of  the  smaller 
domain.  Even  after  the  front  wave  passes  this  boundary,  significant  differences  are 
not  apparent  until  the  first  crest  passes  the  boundary.  The  error  then  starts  to  grow 
and  propagate  back  into  the  domain.  This  differs  from  our  attempts  to  solve  the 
identical  problem  using  the  direct  method — extraneous  wave  reflections  were  noticed 
at  the  boundary  almost  immediately.  This  is  remarkable,  especially  considering 
that  the  non-singular  indirect  method  requires  no  special  treatment  such  ■  one- 
side  derivatives  along  the  edge  of  the  computational  domain.  While  in  the  direct 
method,  numerical  spatial  derivatives,  especially  one-side  derivatives  at  the  boundary 
are  necessary  to  calculate  the  velocities  of  the  nodes.  It  is  known  that  numerical 
differentiation  is  very  likely  to  cause  numerical  instabilities. 

3.4  Conclusions 

The  non-singular  formulation  significantly  reduces  the  time  required  to  compute 
the  influence  matrix  of  the  resulting  algebraic  system.  The  algebraic  systems  are 
still  adequately  well-conditioned  to  allow  effi<aent  iterative  solutions.  Accurate  solu¬ 
tions  can  be  obtained  by  the  desingularized  boundary  method  for  a  large  range  of 
desingulcirization  distances  on  the  order  of  the  mesh  size  (near  =  1).  The  indirect 
method  is  more  efficient  than  the  direct  method.  It  is  easy  to  code  and  takes  less 
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computation  time.  In  addition,  the  indirect  method  performs  better  for  problems 
with  infinitive  boundaries  and  handles  irregular  meshes  more  easily.  Nonlinear  wave 
calculations  using  the  time-stepping  procedure  with  desingularized  boundary  inte¬ 
gral  method  show  promising.  The  method  will  be  examined  in  the  studies  on  more 
complicate  nonlinear  wave  problems  in  the  following  chapters. 
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Figure  3.14;  Waves 
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CHAPTER  IV 


TWO-DIMENSIONAL  SOLITONS  IN  A 
SHALLOW  WATER 


A  disturbance  moving  steadily  with  a  transcritical  velocity  in  a  shallow  water 
can  generate  a  succession  of  solitary  waves,  advancing  upstream  of  the  disturbance. 
A  train  of  weakly  nonlinear  and  weakly  dispersive  waves  develops  downstream  of  a 
region  of  depressed  water  surface  trailing  just  behind  the  disturbance.  Although  the 
disturbance  is  steady  (moving  at  a  constant  speed),  the  flow  generated  is  never  steady 
because  the  flow  is  in  resonance  when  the  disturbance  moves  near  the  critical  speed. 
This  phenomenon  was  first  found  by  Wu  &  Wu  (1982)  using  a  generalized  Boussinesq 
model  for  two-dimensional  long  waves  generated  by  a  moving  surface  pressure  or  bot¬ 
tom  topography  distribution.  This  phenomenon  was  later  experimentally  confirmed 
by  Ertekin,  Webster  &:  Wehausen  (1986)  and  Lee,  Yates  &  Wu  (1989).  It  has  been 
found  that  it  is  important  to  include  the  nonlinear  term  in  the  free  surface  boundary 
conditions  in  order  to  get  the  solitary  waves. 

The  generalized  Boussinesq  model  is  only  good  for  long  waves  (t.e.  e  =  hoj\  is 
small,  where  ho  is  the  mean  depth  of  the  water  and  A  is  a  typical  wave  length).  A  high- 
order  equation  is  given  in  Dommermuth  &  Yue  (1987)  for  the  problem.  For  a  weak 
disturbance  near  the  critical  speed,  a  simplified  forced  KdV  equation  (fKdV),  derived 
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from  the  Boussinesq  equation,  can  be  used  to  model  the  phenomenon  (Akylas  1984, 
Cole  1985,  Wu  1987,  and  Lee  c/  al.  1989).  Attempts  have  been  made  to  extend  the 
formulations  to  three-dimensional  waves  with  side  walls  (Mei  1986,  Katsis  &  Akylas 
1987,  Ertekin  et  al.  1986,  and  Pederson  1988). 

Most  theoretical  studies  on  this  phenomenon  are  based  on  perturbation  methods. 
Recently,  Choi,  tt  al.  (1990)  used  a  finite  element  method  based  on  Luke’s  variational 
principle  for  the  generation  of  solitons  by  a  ship  moving  in  a  channel.  Grilli  & 
Svendsen  (1989)  used  the  mixed  Eulerian-Lagrzmgian  time  stepping  formulation  with 
a  boundary  integral  equation  technique  (singulair  formulation)  for  a  related  problem 
studying  the  runup  and  reflection  of  a  two-dimensional  solitary  wave  in  a  numerical 
tank.  In  their  method,  GriUi  Svendsen  solve  the  Laplace  equation  numerically 
with  the  fully  nonlinear  free  surface  boundary  conditions.  These  numerical  methods 
have  broader  applications  than  the  perturbation  methods. 

In  this  chapter,  the  waves  generated  by  a  free  surface  pressure  or  a  bottom  to¬ 
pography  moving  at  the  critical  speed  are  studied  using  both  the  fKdV  model  (Wu, 
1987)  and  the  fully  nonlinear  model  with  the  desingularized  method.  The  results 
from  the  two  different  models  are  compared.  The  difference  in  the  waves  due  to 
the  free  surface  pressure  and  the  bottom  topography  with  an  identical  dimensionless 
distribution  is  studied  by  the  fully  nonlinear  model.  In  the  fKdV  model,  these  two 
different  disturbances  will  generate  the  same  waves  (Wu  1987).  We  will  also  study 
the  nonlinear  waves  due  to  submerged  bodies,  which  are  difficult  for  the  perturbation 
methods  mentioned  above  to  consider.  The  bottom  is  horizontally  flat  except  in  the 
area  of  the  varying  bottom  topography.  The  free  surface  extends  to  infinity  at  both 
ends.  Fig.  4.1.  The  flow  is  assumed  irrotational.  The  problem  is  nondimensionlized 
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Figure  4.1:  Problem  definition  for  solitons  in  shallow  water 


by  the  mean  water  depth  ho-,  the  gravity  acceleration  g  and  the  fluid  density  p. 

4.1  fKdV  Equation 

The  fKdV  equation  for  a  free  surface  pressure  and  a  bottom  topography  distri¬ 
bution  can  be  derived  from  the  Boussinesq  equations  (Wu  1987)  and  has  the  form 

3  11  1 

“  (1  d"  2^)Vx  ~  ~  2^^**  ~  2^ (^•^) 

where  pa  =  Pa(x  -f  F„t)  and  b  =  6(z  -I-  F„t)  represent  a  free  surface  pressure  and 
a  bottom  topography  moving  in  — z  direction,  respectively,  and  F„  is  the  Froude 
number  based  on  the  mean  water  depth.  Given  the  forcing  F(z  -f-  F„t)  =  pa(z  -I- 
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Fnt)  +  b{x  +  Fnt)  associated  with  the  initial  conditions,  Eq.  (4,1)  can  be  solved  for 
the  wave  elevation  rj.  Detail  of  the  derivation  can  be  found  in  Wu  (1987). 

Within  the  fKdV  model,  the  free  surface  pressure  and  the  bottom  topography 
disturbance  play  the  same  forcing  role,  i.e.  the  waves  generated  by  a  pressure  dis¬ 
turbance  will  be  the  same  as  the  waves  due  to  a  topography  disturbance  if  the  two 
disturbances  have  the  same  distribution  and  magnitude.  This  implies  that  in  an 
experiment,  a  free  surface  pressure  distribution  can  be  replaced  by  the  equivalent 
topography  bottom  that  is  easier  to  control  (Lee  et  al.  1989).  However,  they  are  not 
interchangeable  if  the  parameter  ranges  are  beyond  the  fKdV  assumptions. 

No  closed  form  solution  to  (4.1)  has  been  found.  A  numerical  approach  is  nec¬ 
essary.  For  the  concern  of  numerical  stability  in  solving  the  fKdV  equation,  it  is 
preferable  to  replace  the  dispersion  term  i/x,*  in  (4.1)  by  based  on  the  lowest 
order  solution.  Thus  (4.1)  becomes, 

3  11 

-  (1  +  ^  2^*' 

This  was  first  pointed  out  by  Benjamin  et  al.  (1972)  and  later  used  by  Wu  (1987) 
and  Lee  et  al.  (1989) 

Following  Wu  (1987),  we  solve  the  fKdV  equation  (4.2)  numerically  with  a  time¬ 
stepping  procedure  based  on  the  Euler’s  predictor-corrector  algorithm  to  advance 
time  without  iteration  at  the  correction  stage.  A  second-order  central-difference  ap¬ 
proximation  is  used  for  the  spatial  derivatives.  Since  this  is  an  implicit  method,  both 
the  prediction  and  correction  stages  reduce  to  the  inversion  of  a  constant  tridiagonal 
matrix,  which  is  strictly  diagonally  dominant. 

A  computational  window  moving  with  the  disturbance  is  employed.  The  time 
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step  At  and  the  shifted  distance  Ax  (usuedly  the  mesh  size)  are  chosen  such  that 

Ax 


J  = 


F^At 


(4.3) 


is  an  integer  ( J  =  1, 2, 3, ...).  Thus  the  window  is  shifted  every  J  time  steps  and  if  we 
view  the  solution  only  when  the  window  is  shifted,  the  position  of  the  disturbance  in 
the  window  does  not  change.  Each  time  the  window  is  shifted  one  point  is  dropped 
at  the  downstream  side  of  the  window  and  a  new  point  emerges  into  the  window 
from  upstream.  Because  of  the  finite  computational  domain,  a  treatment  for  the 
open  boundairy  conditions  is  necessary  to  reduce  non-physical  reflections  from  the 
the  boundaries.  Wu  &  Wu  (1982),  Wu  (1987)  and  Lee  et  al.  (1989)  use 


m  =  ±Vt  (4.4) 

where  the  “-f-”  is  taken  for  the  upstream  while  ”  for  the  downstream.  Although 
this  open  boundary  condition  is  only  appropriate  for  lineair  waves  with  unit  phase 
speed,  it  has  been  shown  that  it  can  eliminate  most  of  the  non-physical  small  am¬ 
plitude  wave  reflection  from  open  boundairies.  The  open  boundary  condition  (4.4)  is 
discretized  by  a  forward  difference  in  time.  For  the  spatial  discretization,  a  backward 
difference  is  used  for  the  upstream  point  and  a  forward  difference  for  the  downstreaun 
point.  Even  though  the  open  boundary  condition  (4.4)  is  used,  the  value  of  rj  for  the 
new  coming  point  upstream  still  needs  to  be  provided.  It  is  not  clearly  stated  how 
the  value  of  tj  for  the  coming  point  is  provided  in  the  above  mentioned  papers.  We 
provide  the  value  by  quadratic  extrapolation  using  the  three  closest  points. 


4.2  Fully  Nonlinear  Model 

We  apply  the  indirect  desingular ized  method  to  this  problem.  A  submerged 
moving  body  can  also  be  considered  in  this  method.  At  each  time  instant,  the 
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potential  is  constructed  as, 

Nt  JV* 

y)  =  H  y,  x/,  y/)  +  53  y,  Xk'',yk'"),  (4.5) 

i=i  fc=i 

where  (x,  y)  is  the  field  point.  is  the  potential  due  to  the  source  at  point  (x/,y/) 
above  the  free  surface  and  its  image  source  about  the  z  =  —1  plane,  and  the 
potential  due  to  the  source  at  point  inside  the  rigid  surface  and  its  image, 

=  log[(x  -  x/)*  +  (y  -  y/f]  +  log[(x  -  x/f  +  (y  -  y/  +  2)^],  (4.6) 

=  log[(x  -  x*'‘)^  +  (y  -  y*'^)^]  +  log[(x  -  x**)^  +  (y  -  yt^  +  2)^].  (4.7) 

The  function  4>{x,y)  automatically  satisfies  the  Laplace  equation  aind  the  boimdary 
condition  on  the  bottom  except  the  area  of  a  bottom  topography.  The  Dirichlet 
condition  on  the  free  surface  and  the  Neumann  condition  on  other  rigid  body  surfaces 
remain  to  be  satisfied.  The  source  strength  and  o-**  are  determined  by  satisfying 
the  boundary  conditions  at  the  chosen  control  points  on  the  free  smrfaice  and  rigid 
surface, 

si  s'' 

'^<r/^^{xi,yi,x/,y/)  +  53  <T\$'^(xi,yi,Xfc\yjfe'‘)  =  ^o(i.,y,) 
j=i  fe=i 

(i  =  l,2,...,yv0,  (4.8) 

si  a  a 

'^<r^j-^^Hxhyhx/,y/)  +  '^<T^k-^^’*{xi,yi,Xk^,yk'')  =  Vh  •  n(x/,y,) 

j=i  *=i 

(/=l,2,...,iV*),  (4.9) 

where  (x,-,  y,)  is  the  control  point  on  the  free  surface  and  (x/,  yi)  is  the  control  point  on 
the  rigid  surface.  <j>o  is  known  from  the  previous  time.  Equations  (4.8)  cind  (4.9)  form 
a  system  of  {N^  +  N^)  x  {N^  +  N*')  linear  equations.  After  solving  (4.8)  and  (4.9) 
for  and  <Tk^,  the  flow  is  determined  and  the  flow  velocities  can  be  calculated  by 
taking  the  spatial  derivative  of  (4.5).  The  free  surface  position  and  the  flow  potential 
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on  it  can  then  be  updated.  The  pressure  on  the  body  can  also  be  calculated  using 
(2.16)  and  integrating  the  pressure  gives  the  forces  on  the  body. 

The  fourth-order  Runge-Kutta-Fehlberg  method  and  the  same  computational 
window  shifting  technique  for  the  fKdV  equation  is  used  in  the  time  stepping.  Each 
time  the  window  is  shifted  one  node  point  is  dropped  at  the  downstre<im  side  of 
the  window  without  any  treatment  since  no  spatial  derivative  of  the  &ee  surface 
is  required  in  this  algorithm.  For  the  upstream  boundary,  the  position  and  the 
potential  of  the  incoming  point  are  pre-computed  using  some  extra  upstream  points 
outside  the  window,  which  are  convected  by  the  computed  velocity  from  the  already 
solved  source  strength.  The  extra  points  do  not  contribute  to  the  discretized  integral 
equation. 

4.3  Numerical  Results 

For  the  fully  nonlinear  method,  the  effects  of  the  mesh  size,  the  error  telorance 
for  the  GMRES  and  the  error  tolerance  for  the  Runge-Kutta-Fehlberg  subroutine 
on  the  solution  have  been  numerically  tested.  It  has  been  found  that  as  long  as  the 
GMRES  error  tolerance  is  less  than  +  N^),  the  Runge-Kutta-Fehlberg  error 

telorance  less  than  10“®,  zmd  if  the  mesh  has  at  least  15  nodes  per  wave  length,  the 
solution  will  be  adequately  resolved. 
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4.3.1  Solitons  due  to  a  free  surface  pressure 


The  waves  generated  by  a  moving  surface  pressure  distribution  are  studied.  The 
pressure  distribution  (Wu  1987)  is  used, 

3Pm{H-cos[^(x  +  F„<)]}  |x  +  F„t|<f 
0  |x  +  Fnt\  >  f 

This  will  also  check  our  programs  for  the  fKdV  model  and  the  fully  nonlinear  model. 


P(x  +  Fnt)  =  p„(x  +  Fnt)  =  < 


First  the  results  by  the  fKdV  and  the  fully  nonlinear  algorithm  for  a  relatively 
weak  forcing  (Pm  =  0.02)  are  compared.  The  results  using  301  node  points  and  a 
time  step  At  =  0.2  in  the  two  models  are  shown  in  Fig.  4.2.  No  upstream  runaway 
soliton  is  observed  since  the  disturbance  is  weak  and  the  period  for  the  generation 
of  the  solitons  is  very  large.  In  general,  the  results  of  the  two  methods  are  similar 
although  the  transient  trailing  waves  computed  by  the  fKdV  program  seem  to  travel 
faster  downstream  away  from  the  disturbance  than  those  by  the  nonlinear  program 
and  the  nonlinear  program  predicts  a  steeper  front  wave.  Both  methods  have  small 
non-physical  downstream  reflections  of  about  the  same  magnitude  after  the  waves 
reach  the  downstream  boundary.  The  results  change  little  with  finer  grids.  It  is 
diflScult  to  isolate  the  causes  of  the  differences  in  the  resiilts  by  the  fKdV  and  the 
nonlinear  computation.  They  may  be  due  to  the  computational  errors  or  to  the 
different  mathematical  models.  It  seems  that  the  difference  is  mainly  due  to  the 
different  models  for  both  the  free  surface  nonlinearity  and  the  two-dimensional  effect. 
This  can  be  seen  by  examining  the  short  time  results  where  the  trailing  waves  have 
not  reached  the  downstream  boundary  and  the  solutions  axe  not  affected  by  the 
non-physical  reflection. 


For  stronger  forcing,  Pm  —  0.1,  both  models  predict  the  upstream  runaway  soli- 
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tons  shown  in  Fig.  4.3.  The  results  are  in  good  agreement  in  general.  However, 
the  nonlinear  model  again  predicts  steeper  waves  with  a  shorter  period  to  generate 
solitons.  This  becomes  more  obvious  when  the  forcing  is  stronger.  Fig.  4.4  shows  the 
results  for  Pm  =  0.15.  As  seen,  the  nonlinear  model  predicts  a  much  steeper  front 
wave  which  starts  to  break  at  t  =  24.2  and  the  program  stops.  For  Pm  =  0.2,  which 
is  used  in  Wu  (1987)  and  Lee  ei  al.  (1989),  the  waves  by  the  nonlinear  program 
break  sooner  (not  shown  here),  while  the  waves  by  the  fKdV  model  do  not  break. 
This  does  not  necessarily  indicate  that  the  fKdV  model  gives  a  better  solution.  As 
assumed  for  the  derivation  of  the  fKdV  equation  (Wu  1987),  the  forcing  Pm  must 
be  O(e^).  A  typical  wavelength  can  be  estimated  for  the  Ccise  Pm  =  0.1  shown  in 
Fig.4.3  from  the  figure.  The  wavelength  A  is  foimd  approximately  10.  Therefore, 
eisapproximatelyO.l.  So  the  forcing  Pm  should  be  0(10—*).  The  values  for  Pm  used 
in  Wu  (1987)  seem  far  beyond  this  limitation.  Wave  breaking  was  observed  in  the 
experiments  (Lee  et  al.  1989). 

Fig.  4.5  shows  the  results  for  a  surface  pressure  at  the  supercritical  speed  (F„  = 
1.5).  Both  models  reveal  the  similar  characters  of  the  flow  motion.  The  dowmstream 
trcinsient  waves  travel  at  slower  speeds  than  the  disturbance  and  eventually  only  a 
single  steady  soliton  remains  traveling  at  the  same  speed  as  the  forcing.  Also,  since 
the  flow  is  not  in  resonance,  the  soliton  height  is  smaller  than  those  generated  by 
the  disturbances  near  the  critical  speed. 

4.3.2  Comparison  of  free  surface  pressure  and  bottom  topography 

In  the  fKdV  model,  the  waves  due  to  a  pressure  disturbance  are  the  same  as 
the  waves  due  to  a  topography  disturbance  if  their  dimensionless  distributions  are 
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identical.  This  is  due  to  the  fact  that  the  £KdV  model  assumes  small  c  and  uses 
the  vertical  average  of  the  flow  velocity  over  the  water  layer,  thus  not  taking  into 
account  the  two-dimensional  effects  of  the  vertical  variation.  The  differences  in  the 
waves  generated  by  the  free  surface  pressure  and  the  bottom  topography  of  identical 
distribution  using  the  fully  nonlinear  model  are  examined. 

Fig.  4.6  shows  the  waves  due  to  a  moving  semi-elliptical  bottom  topography  with 
the  major  axis  L  and  minor  axis  bo  being  2.0  and  0.2,  respectively.  Fig.  4.7  shows 
the  waves  generated  by  the  free  surface  pressure  of  the  same  dimensionless  elliptical 
distribution.  The  difference  in  the  waves  due  to  the  two  different  disturbances  are 
very  significant.  The  free  surface  pressure  behaves  as  a  stronger  forcing  than  the 
bottom  topography.  The  waves  due  to  the  pressure  start  to  break  at  about  t  =  55. 
The  breaking  occurs  at  the  first  wave  peak  downstreaun  of  the  disturbance. 

4.3.3  Solitons  due  to  submerged  cylinders 

Waves  generated  by  submerged  cylinders,  which  cannot  be  modeled  by  the  fKdV 
equation  but  can  be  easily  handled  by  the  nonlinear  model,  are  calculated.  Fig.  4.8 
shows  the  waves  due  to  a  circular  cylinder  with  the  radius  R  =  0.15  and  the  depth 
of  the  centroid  hg  =  0.6  moving  at  F„  =  1.0.  Solitons  generated  are  clearly  seen. 
Fig.  4.9  shows  the  waves  generated  by  a  larger  circular  cylinder  with  R  =  0.2  at 
the  same  speed.  The  the  cylinder  with  larger  radius  implies  a  stronger  disturbance 
and  therefore  produces  larger  waves  and  the  first  wave  downstream  starts  to  break 
at  about  t  =  45.  Fig.  4.10  shows  the  waves  due  to  the  cylinder  {R  =  0.15)  at 
the  supercritical  speed  (F„  =  1.5).  Again  only  one  single  soliton  eventually  moves 
with  the  disturbance  as  in  the  case  of  the  surface  pressure  (Fig.  4.5).  Fig.  4.11  and 
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Fig.  4.12  show  the  waves  due  to  elliptical  cylinders.  In  the  former  case,  the  ratio  of 
the  major  aocis  to  the  minor  axis  is  RyfRx  =  0.075/1.0.  In  the  later  case,  the  cyhnder 
has  a  larger  minor  axis  {Ry/Rx  =  0.1/1. 0)  representing  a  stronger  disturbance.  The 
breaking  of  the  downstream  wave  occurs  in  the  later  case. 

4.4  Conclusions 

The  fully  nonlinear  model  predicts  the  runaway  solitons  for  three  types  of  dis¬ 
turbances  moving  near  the  critical  speed:  surface  pressure,  bottom  topography  and 
submerged  cylinder.  The  numerical  results  by  the  fully  nonlinear  model  agree  with 
the  fKdV  results  for  the  weak  surface  pressure  disturbance.  For  the  strong  dis¬ 
turbance,  the  fuUy  nonlinear  model  predicts  larger  solitons  than  the  fKdV.  These 
sohtons  are  so  steep  that  they  break  either  upstream  or  downstream.  According  to 
the  fully  nonlinear  calculations,  for  an  identical  distribution,  the  free  surface  pressure 
generates  larger  waves  than  the  bottom  topography  and  the  difference  is  significant. 


54 


y _ P(a:.  t) 


Pm  =  0.02 

^  ^  ^  r  *  f  *  r-r**f  r-r  _  2.0 


Figure  4.2:  Waves  due  to  free  surface  pressure  (P^  =  0.02,  F„  =  1.0,  At  =  0.2) 


Figure  4.3:  Waves  due  to  free  surface  pressure  =  0.1,  =  1.0,  At  =  0.2) 
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Figure  4.4:  Waves  due  to  free  surface  pressure  {Pm  =  0.15,  Fn  =  1.0,  Af  =  0.2) 
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Figure  4.5:  Waves  due  to  free  surface  pressure  =  0.1,  F„  =  1.5,  At  =  0.2) 
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Figure  4.9:  Waves  due  to  submerged 
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CHAPTER  V 


WAVES  GENERATED  BY  A  SUBMERGED 

SPHEROID 


In  this  chapter,  the  desingularized  boundary  integral  equation  is  applied  to  a 
spheroid  moving  below  a  free  surface.  The  waves  due  to  a  submerged  spheroid  have 
been  studied  by  many  researchers  and  used  to  check  various  theories  and  numerical 
methods  (Doctors  &  Beck  1987,  Havelock  1931,  Farell  1973a,b,  Bertr<im  1990,  and 
Magee  1991). 

For  a  steady  problem,  Doctors  h  Beck  (1987)  lineau-ize  the  free  surface  condition 
about  uniform  flow  and  solve  the  resulting  classical  Neumann- Kelvin  problem  using 
a  panel  method  with  Kelvin  sources  of  constant  strength  over  each  panel.  Bertram 
(1990)  uses  the  Rankine  source  method  of  Jensen,  Mi  &  Soding  (1986).  In  their 
method,  discrete  sources  are  distributed  above  a  finite  section  of  the  free  surface  to 
fulfill  the  free  surface  boundary  conditions.  On  the  body  surface,  polygonal  panels 
of  constant  strength  are  distributed.  The  nonlinear  free  surface  boimdary  condition 
is  met  using  an  iterative  scheme  that  linetirizes  around  an  approximation  from  pre¬ 
vious  iteration.  The  potential  and  wave  elevation  are  improved  in  each  iteration. 
The  radiation  and  open  boundary  conditions  are  enforced  by  shifting  sources  versus 
collocation  points  on  the  free  surface. 
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For  an  unsteady  problem  with  forward  speed,  Magee  (1991)  uses  a  time-domain 
analysis  to  solve  the  unsteady  Neumann- Kelvin  problem.  In  his  analysis,  the  body 
boundary  condition  is  applied  on  the  exact  time-dependent  body  surface,  while  the 
linear  free  surface  condition  is  retained  which  allows  the  use  of  the  time-domain 
Green  function.  Since  the  Green  function  satisfies  the  free  surface  and  radiation 
conditions,  unknowns  need  only  be  distributed  on  the  body  surface  at  any  instant  of 
time.  However,  because  of  the  memory  effect  from  the  free  surface,  a  large  computer 
storage  is  required  to  store  the  strengths  and  the  positions  of  the  sources  in  order 
to  evaluate  the  convolution  integral.  The  storage  requirement  is  proportional  to  the 
simulation  time,  which  may  prohibit  large  time  simulations. 

We  apply  the  desingulaxized  method  for  the  unsteady  problem.  Discrete  Rankine 
sources  are  distributed  above  the  truncated  free  surface  and  inside  the  spheroid. 
The  evaluation  of  Rankine  sources  is  very  fast,  especially  on  a  supercomputer  since 
a  very  high  level  of  vectorization  can  be  achieved  because  of  the  simplicity  of  the 
Rankine  source.  The  errors  due  to  the  free  surface  truncation  is  expected.  With  a 
proper  implementation  of  the  open  boundary  conditions,  the  truncation  error  can  be 
reduced  to  minimal. 

The  spheroid,  with  the  diameter-length  ratio  Dj L  equal  to  0.2,  starts  to  move 
horizontally  from  rest  in  the  -i  direction  at  the  depth  H  below  the  free  surface  and 
reaches  an  ultinnate  steady  velocity.  The  flow  is  at  rest  until  the  spheroid  starts  to 
move  and  will  eventually  reach  a  steady  state  ir.  the  moving  coordinate  system.  The 
flow  is  sj  mmetric  about  the  Oxz  plane.  The  computation  will  take  advantage  of  this 
symmetry.  The  problem  is  nondimensionized  by  making  the  length  of  the  spheroid 
L,  the  gravity  and  the  fluid  density  p  unity. 
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Figure  5.1:  The  coordinate  and  the  discretization  of  the  spheroid 

5.1  Numerical  Aspects 

The  free  surface  is  discretized  using  Nl  nodes  and  nodes  in  the  x  «ind  y  direc¬ 
tions,  respectively,  to  form  =  iV/  x  free  surface  nodes  within  a  computational 
domain  0  <  x  <  72,  and  0  <  y  <  Ry.  The  initied  grid  has  an  equal  spacing  in  the 
X  direction  and  an  algebraically  increasing  spacing  in  the  y  direction  further  from 
the  symmetric  plane.  On  the  spheroid,  the  grid  lines  are  spaced  uniformly  in  the 
circumferential  direction.  The  grid  lines  have  a  cosine  spacing  in  the  longitudinal 
direction.  The  body  has  a  grid  with  elements  in  x  direction  and  elements 
in  the  circumferential  direction,  resulting  in  TV*  =  (TV*  —  1)  x  (TV|  -f  1)  -b  2  nodes 
including  the  two  end  points  on  the  body  surface,  Fig.  5.1. 

We  constrttct  the  solution  using  a  source  distribution  on  a  surface  {S/)  above 
the  free  surface  and  a  source  distribution  on  a  surface  inside  the  body  surface. 
For  a  better  computational  open  boundary  behavior,  we  add  the  negative  images  of 
the  source  distribution  on  5*  about  the  undisturbed  free  surface.  So  the  potential  is 
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expressed  as 


<t>{x)  = 


I  /v  + 1 L  -  ]i:bi 


I  -  X'i 


X  —  X.- 


where  is  the  point  on  5/',  x*  is  the  point  on  Sh  and  I*  is  the  image  of  x*. 
By  applying  the  free  surface  boundary  condition  (2.12)  and  the  body  boundary 
condition  (2.13),  we  obtain  a  system  of  linear  equations  for  the  unknown  strength  of 
the  singularities, 

N>  1  Af* 

— ;7t)  + 

J  =  1 

I  •  -  I  I  -  -  I 

(i  =  l,2,...,iV^)  (5.2) 


1 


1  ,  ^J\i  '  ^''1  l|  ^  -6!  I  -*  ~ 


Nt  a 


1 


Xii  -  x< 


1 


1 


arfcf-x^l  |x6,-x>| 

(i  =  l,2,...,JV^) 


(5.3) 


where  x”/-  and  xjj  are  the  control  points  on  5/  amd  S/,.  After  cr^  and  are  determined, 
the  fluid  ptirticle  velocities  on  the  free  surface  can  be  calculated.  Then  the  time¬ 
stepping  procedure  integrates  the  free  surface  boundary  conditions. 


Similar  to  the  two-dimensional  problem  in  Chapter  4,  the  fourth-order  Runge- 
Kutta-Felhberg  method  and  the  computational  window  shifting  technique  are  used 
in  the  time  stepping  with  a  simpler  implementation  of  the  open  boundary  conditions. 
Each  time  the  window  is  shifted  one  row  of  node  points  are  dropped  at  the  down¬ 
stream  side  of  the  window  without  any  special  treatment  since  no  spatial  derivative 
of  the  free  surface  is  required.  For  the  upstream  boundary,  a  row  of  new  node  points 
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enter  the  window  with  homogeneous  conditions.  Because  of  an  additional  dimension, 
the  flow  decays  from  the  disturbance  more  rapidly  in  the  three-dimensional  case  than 
in  the  two-dimensional  case.  Therefore,  the  easier  and  simpler  implemention  of  the 
open  boundary  condition  can  be  used  in  mainy  three-dimensional  wave  problems.  It 
will  be  shown  that  this  technique  works  very  well  in  the  numerical  computations. 

The  pressure  on  the  body  surface  is  evaluated  using  the  Bernoulli  equation: 

-P  =  ^  +  +  ^  =  ^  +  (5  V  ^  -  K)  •  (5-4) 

where  ^  =  (^  +  K  is  the  substantial  derivative  of  the  potenticil  at  fixed  points 
on  the  body  surface.  The  nondimensional  velocity  of  the  spheroid  is 

Vi  =  -Fnil  -  e-'")  I  (5.5) 

The  Froude  number  is  defined  as  Fn  =  U where  U  is  the  ultimate  steady 
velocity  of  the  spheroid. 


The  drag  hft  and  moment  M  on  the  body  are  calculated  by  integrating  the 
pressure  over  the  body  surface  using  (2.17)  and  (2.18)  in  conjunction  with  Simpson’s 
rule.  The  hydrodynamic  coefficients  are  defined  as 


Cw  =  — 


C.  =  T 


F, 


\pU'^A  '  ^  \pU^A  ' 

where  A  is  the  surface  area  of  the  spheroid. 


and  C\f  = 


M 


\pmAL  ’ 


(5.6) 


The  mesh  size  on  the  free  surface  is  usually  larger  them  that  on  the  body  by  about 
10  times  and  the  differences  in  the  magnitudes  of  the  coefficient  matrix  elements  are 
large.  Therefore,  the  resulting  system  for  cr/  and  <T(,,  (5.2)  and  (5.3),  as  a  whole 
is  likely  to  be  poorly  conditioned.  To  avoid  this  numerical  difficulty,  we  split  the 
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system  into  two  subsystems,  (5.2)  for  <7/  and  (5.3)  for  which  are  alternately 
solved  using  LU  decomposition  for  each  subsystem.  In  other  words,  <7/  is  solved  by 
(5.2)  with  <76  known  and  <7^  is  then  corrected  by  (5.3)  with  <Tj  known.  The  process 
goes  forth  and  back  between  (5.2)  and  (5.3)  until  the  error  criterion  is  met.  Each  set 
of  equations  is  much  better  behaved  and  more  accurate  solutions  can  be  obtained. 
Another  advantage  of  splitting  is  that  the  coefiBicient  matrix  for  does  not  change 
with  time  and  needs  only  to  be  inverted  once  for  the  entire  time  simulation.  The 
matrix  for  <7/  does  not  change  during  the  iteration  between  the  body  auad  the  free 
surface  and  only  needs  to  be  inverted  once  for  the  current  instant  of  time.  Of  course, 
the  coefficient  matrix  for  <7/  changes  at  the  next  instant  of  time. 

The  pressure  on  the  body  is  evaluated  at  the  node  points.  The  substantial  deriva¬ 
tive  of  the  potential  ^  in  (5.4)  is  calculated  using  a  four-point  forward  difference 
scheme.  The  fourth-order  Runge-Kutta-Fehlberg  method  is  used  in  the  nonlinear 
free-stirface  integration.  An  initial  time  increment  is  set,  but  is  modified  by  the 
Runge-Kutta-Fehlberg  subroutine  where  appropriate. 

5.2  Results 

On  the  free  surface,  we  use  61  x  16  nodes  within  0  <  x  <  7.5  and  0  <  y  <  1.875, 
and  symmetry  about  y  =  0  is  assumed.  The  spacing  in  the  y  direction  increases  by 
10  percent  for  each  row  of  nodes  further  from  the  centerline.  For  comparisons  to  the 
results  presented  in  Doctors  k  Beck  (1987),  we  use  the  same  submergence  depths  of 
the  spheroid  {H f  L  =  0.16  <ind  0.245). 

Fig.  5.2  shows  isometric  and  contour  plots  of  the  waves  caused  by  the  spheroid 
for  H/L  =  0.245  and  Fr  =  0.6  at  t  =  25.  A  smooth  startup  /x  =  2  in  (5.5)  and 
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NI  =  2Ng  =  16  axe  used.  Fig.  5.3  shows  the  drag,  lift  and  moment  acting  on  the 
spheroid  as  a  function  of  time.  They  approach  the  steady  state  values  after  the 
body  has  moved  10  to  15  body  lengths.  Fig.5.4  shows  the  effects  of  two  different 
startups  of  the  body  on  the  hydrodynamic  forces.  In  the  two  startups  (/z  =  2  and 
oo),  the  body  experiences  large  differences  in  force  acting  on  it  during  the  trainsition; 
eventually  the  forces  approach  the  same  steady  state  values.  We  also  notice  that  the 
body  experiences  a  negative  drag  for  a  short  time  soon  adter  the  impulsive  startup. 

Fig.  5.5  compares  the  pressure  on  the  spheroid  using  the  present  method  and  the 
Neumann-Kelvin  calculation  (Magee  1990)  at  t  =  25  for  the  same  conditions.  The 
compau'ison  is  made  for  the  pressure  on  the  top,  bottom  and  side  of  the  spheroid.  In 
the  Neumann-Kelvin  calculation,  the  body  surface  is  divided  into  flat  panels  with 
distributions  of  constant  source  strength.  The  strength  is  determined  by  satisfying 
the  body  boundary  condition  at  the  centers  of  the  panels.  The  pressure  is  calculated 
at  the  center  points.  The  pressure  elsewhere  is  obtained  by  interpolation.  The 
differences  between  the  linear  cind  the  nonlinear  results  are  noticeable. 

Fig.  5.6  shews  the  convergence  of  the  drag,  lift  and  moment  with  respect  to  the 
node  number  with  =  27V|  =  8,  12,  16  and  20.  For  all  these  cases,  the  free 
surface  grid  (61  x  16)  adequately  resolves  the  waves. 

Comparisons  of  the  hydrodynamic  coefficients  at  different  Froude  numbers  to 
the  Neumman- Kelvin  calculations  (Doctors  &  Beck  1987)  and  other  calculations  by 
Bertram  et  al  (1991)  and  Boppe  (1990)  are  shown  in  Fig.  5.7  for  Hj L  =  0.245  (Note: 
the  moments  in  the  Fig.  5.7  are  taken  about  the  vertical  projection  of  the  spheroid 
centroid  onto  the  undisturbed  free  surface.  The  moments  shown  in  the  other  figures 
are  taken  about  the  centroid  of  the  spheroid).  Bertram’s  calculation  is  based  on  the 


72 


method  of  Jensen,  Mi  and  Soding  (1986).  Boppe’s  calculation  uses  a  Dawson-type 
scheme.  Their  calculations  solve  the  boundary  value  problem  for  steady  flow.  Our 
cadculations  solve  the  initial  boundary  value  problem  for  unsteaidy  flow  generated  by 
the  spheroid  starting  from  rest.  After  the  spheroid  reaches  it  final  constant  speed,  the 
flow  will  become  steady  as  <  — ♦  oo.  Our  results  in  Fig.  5.7  are  obtained  by  curve  fit 
assuming  the  unsteady  components  decay  at  0{t~^)  according  to  Wehausen  (1964). 
A  finer  free  surface  grid  (71  x  16)  is  used  for  smaller  Froude  numbers  (F„  =  0.3  and 
0.4)  to  resolve  the  shorter  waves.  Our  results  agree  well  with  the  linear  results  of 
Doctors  &  Beck  (1987)  at  low  Froude  number  bs  expected  since  the  nonlinear  effects 
are  presumably  small  for  low  Froude  number.  As  the  Froude  number  increases, 
the  difference  between  our  results  and  the  linear  results  increase  since  the  nonlinecir 
effects  are  expected  to  become  larger.  Also  notice  that  there  are  differences  between 
different  calculations,  especially  for  large  Froude  number.  This  may  be  due  to  the 
different  truncations  in  our  calculations  and  the  other  calculations  or  due  to  the 
error  from  the  curve  fitting  for  large  time  asymptotic  values  in  our  calculations. 
Nevertheless,  the  agreement  is  good  in  general. 

For  HjL  =  0.16,  nonlinear  effects  seem  very  strong.  For  all  attempted  Froude 
numbers  (i^„=0.4,  0.5,  0.6,  0.7  and  0.8),  the  free  surface  in  our  computation  is 
sucked  down  and  touches  the  body  surface  which  in  turn  stops  the  computation.  A 
smoother  start  up  could  delay  the  surface  piercing  of  the  body.  It  seems  that  in 
our  calculations,  the  acceleration  of  the  body  specified  is  too  large  so  that  the  initial 
start-up  causes  surface  piv,fcing  of  the  body.  It  is  also  possible  for  this  to  happen 
even  if  the  acceleration  is  very  small  but  the  final  speed  is  high.  Surface  piercing  of 
the  body  also  occurs  in  Bertram’s  steady  calculation  (Bertram  1990)  at  high  Froude 
number  {Fn  =  0.8).  The  linear  calculations  arc  not  affected  by  this  because  the  free 
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surface  boundary  conditions  are  satisfied  at  2  =  0.  The  surface  piercing  of  the  body 
needs  to  be  verified  by  experiments. 

The  waves  generated  by  the  spheroid  and  those  generated  by  the  relevant  simple 
source-sink  pair  disturbance  are  compared  (Fig.  5.2  and  Fig.  5.8).  The  strength 
of  the  pair  and  the  distance  between  the  source  and  sink  are  determined  to  give  a 
Rankine  ov2tl  having  the  same  length  and  midsectional  area  as  the  spheroid  moving 
in  an  infinite  fluid  The  comparison  becomes  meaningless  if  the  disturbance  is  too 
close  to  the  free  surface  because  the  simple  source-sink  no  longer  represents  the 
body  well.  The  same  depth  of  submergence,  location  of  the  center  and  motion  of  the 
disturbance  as  those  for  the  spheroid  are  used  for  a  direct  comparison.  Compcirison 
of  Fig.  5.2  and  Fig.  5.8  shows  that  the  wave  patterns  of  the  spheroid  and  the  relevant 
source-sink  are  very  similar  except  that  the  spheroid  generates  steeper  waves  near 
the  stern. 

The  computations  are  carried  out  on  a  Cray  Y-MP.  The  results  shown  in  Fig.  5.2 
take  approximately  4.3  CPU  seconds  to  solve  the  boundary  value  problem  (5.2) 
and  (5.3)  in  which  10  percent  of  this  time  is  required  for  the  matrix  setup.  About 
5  iterations  between  (5.2)  and  (5.3)  are  required.  The  CPU  time  for  the  entire 
time  simulation  is  45  minutes.  For  the  simpler  disturbance  in  Fig.  5.8  it  takes 
approximately  3.9  CPU  seconds  to  solve  (5.2)  and  40  minutes  for  the  entire  time 
simulation.  In  the  later  caise,  only  (5.2)  needs  to  be  solved.  The  difference  in  the 
CPU  time  for  the  two  computations  is  not  significant  because  the  splitting  procedure 
for  the  body  problem  only  requires  an  additional  evaluation  of  multiple  right-hand 
sides. 


A  faster  iterative  matrix  solver  can  be  used  effectively  for  the  source-sink  pair 
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disturbance  to  further  reduce  the  computation  time.  Although  GMRES  could  have 
been  used  for  the  whole  system  in  the  spheroid  example,  it  would  require  special 
preconditioning.  We  tried  solving  the  whole  system  for  the  spheroid  example  using 
both  LU  decomposition  and  GMRES  without  preconditioning.  The  LU  algorithm 
gave  inaccurate  solutions  while  GMRES  did  not  converge. 

Fig.  5.9  shows  the  effect  of  desingularization  (/^  =  0.5, 1.0, 2.0, 3.0  and  4.0)  on 
the  wave  computations.  The  results  are  not  significantly  affected  except  for  U  =  0-5, 
which  becomes  too  small  and  the  integration  by  summation  (or  one  point  Gauss 
quadrature)  is  too  crude. 

We  also  performed  an  energy  conservation  check  for  a  fixed  control  volume 
bounded  by  the  free  surface,  the  body  surface,  a  horizontal  bottom  and  four  ver¬ 
tical  surfaces  representing  the  truncated  far-field  boundaries.  The  body  starts  to 
move  from  rest  at  the  center  of  the  control  volume.  The  energy  of  the  fluid  in  the 
control  volume  and  the  work  done  by  the  body  aire  shown  in  Fig.  5.10.  The  energy 
and  the  work  balance  each  other  well  until  t  =  6  when  the  waves  and  the  body  reach 
the  truncated  boundary.  Since  the  energy  flux  is  neglected  in  the  energy  calculation, 
agreement  is  not  to  be  expected  after  the  body  or  the  waves  reach  the  boundary. 

5.3  Conclusions 

Desingularization  performs  well  for  fully  nonhnear  free  surface  problems  of  a 
submerged  spheroid.  Desingularization  is  not  numerically  dispersive  (Fig.  5.9)  nor 
dissipative  (Fig.  5.10).  Iteration  between  the  free  surface  and  the  body  surface 
boundary  conditions  is  required  for  the  not  well-conditioned  matrix  equation  without 
preconditioning.  The  waves  produced  by  the  spheroid  and  the  relevant  source-sink 
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disturbance  are  similar  if  they  are  sufficiently  submerged.  This  indicates  that  a 
simple  disturbance  (with  a  much  simpler  iterative  procedure  and  fewer  unknowns) 
can  be  used  if  the  surface  waves  are  the  main  interest.  In  general,  the  comparison  of 
the  hydrodynamic  coefficients  by  our  results  to  those  by  other  calculations  is  good 
for  HjL  =  0.245.  For  HjL  =  0.16,  our  method  predicts  the  surface  piercing  of  the 
body,  which  needs  to  be  verified  by  experiments. 
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Figure  5.2:  Waves  generated  by  the  spheroid,  above:  isometric  view,  below:  wave 
elevation  contours  (0.02  apeurt) 
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Figure  5.3:  Hydrodynamic  coefficients  vs.  time 
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Figure  5.4:  Influence  of  start-up  of  the  spheroid 
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Figure  5.6:  Sensitivity  of  i 
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Figure  5.8:  Waves  generated  by  the  relevant  source-sink  pair,  above:  i 
below:  wave  elevation  contours  (0.02  apart) 
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CHAPTER  VI 


UNSTEADY  WAKE  AND  INNER- ANGLE 

WAVEPACKETS 


Recent  observations  (SAR  images  from  the  space  shuttle  and  aerial  photographs) 
of  the  ship  wakes  sometimes  reveal  a  narrow  V  wave  pattern  in  addition  to  the 
the  classical  Kelvin  wake  (Fu  &  Holt  1982,  Munk  et  al.  1987,  Brown  et  al.  1989, 
Reed  et  al.  1991  and  the  articles  in  JOWIP  and  SARSEX  Special  Issue,  Journal  of 
Geophysical  Research  1988).  According  to  available  observations,  a  ship  wake  can 
consist  of  the  following  three  types  of  structures,  although  this  classification  is  not 
very  precise:  1)  the  usual  Kelvin  wake  component  consisting  of  transverse  waves 
and  diverging  waves,  both  within  the  classical  19.5  degree  cusp  line;  2)  the  narrow 
V-wake  typically  appearing  in  the  SAR  image  as  a  narrow  bright  V  of  half-angle  2 
to  3  degrees  with  a  centerline  dark  area  within  the  V;  and  3)  the  relatively  isolated 
wavepacket  of  about  10  to  11  degrees,  e.g.  as  recently  observed  in  the  wake  of  the 
U.S.  Coast  Guard  cutter  Point  Brower,  (Brown  et  al.  1989).  The  packet  persists 
for  several  kilometers  behind  the  ship.  Here  we  will  call  this  isolated  wavepacket 
inner-angle  wavepacket  to  distinguish  it  from  the  very  narrow  V-wake  of  half-angle 
2  to  3  degrees.  More  detailed  descriptions  about  the  ship  wake  can  be  found  in  the 
papers  mentioned  above. 
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Although  there  are  many  theoretical  uncertainties  as  well  as  voids  in  experimental 
data,  it  is  believed  that  the  very  narrow  V-wakes  are  associated  with  turbulence, 
internal  waves,  wave  breaking,  nonlinear  effects  as  well  as  surfactant  effects  (Reed 
et  al.  1991).  The  generation  mechanism  responsible  for  inner-angle  waves  and  their 
relation  to  the  ship  characteristics  are  also  not  clear.  The  typical  wavelength  of  the 
inner- angle  waves  is  much  larger  than  those  in  the  very  narrow  V-walces  and  has 
the  same  order  ais  that  of  the  Kelvin  components.  This  can  be  easily  seen  in  the 
aerial  photograph  of  the  wake  of  the  Coast  Guard  cutter  Point  Brower  (figure  2  in 
Brown  et  al.  1989).  The  inner- angle  waves  will  be  studied  in  this  chapter.  The  most 
probable  hypotheses  of  the  causes  for  the  inner-angle  waves  are:  1)  interference  due 
to  the  linear  superposition  of  the  wave  fields  generated  by  the  bow  and  the  stern; 
2)  nonlinear  effects;  and  3)  unsteadiness  of  the  waves.  We  will  investigate  all  three 
cases  in  this  chapter. 

Linear  interference  has  been  studied  by  Hall  &  Buchsbaum  (1990)  using  a  sub¬ 
merged  source-sink  pair  as  well  eis  a  line  sotirce  distribution  to  model  the  ship.  The 
linear  free  surface  boundary  condition  was  satisfied.  They  were  able  to  generate  an 
inner-angle  wavepacket  using  this  model.  However,  the  wavepacket  spreads  faster 
than  the  observed  data.  Brown  et  al.  (1989)  found  that  the  decay  rate  of  the  inner- 
angle  wave  amplitude  predicted  by  linear  steady  Kelvin  wake  theory  was  faster  than 
the  observed  decay  rate.  According  to  Hall  &  Buchsbaum  (1990),  Brown  et  al. 
(1989),  as  well  as  Akylas  et  al.  (1989),  nonlinear  effects  play  an  important  role  in 
the  generation  of  the  inner-angle  waves. 

Most  studies  to  date  on  this  inner-angle  phenomenon  focus  on  nonlinear  effects. 
No  comprehensive  nonlinear  models  have  been  developed  to  calculate  the  entire  wave 
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field  of  a  ship,  although  considerable  study  has  been  conducted  on  certain  aspects. 
Most  studies  use  a  perturbation  method  to  account  for  nonlinear  effects  (Akylas 
et  al.  1988,  Brown  et  al.  1989,  and  Hall  &  Buchsbaum  1990).  A  two-dimensional 
nonlinear  Schrodinger  equation  for  the  wave  envelope  is  derived  assuming  a  nsnrow 
band  of  frequencies  for  the  carrier  waves  adong  the  ray  of  the  inner-angle  wavepacket. 
The  wave  envelope  from  the  nonlinear  Schrodinger  equation  can  support  solitary 
wave  solutions  (hence  the  terminology  iimer-angle  soli  ton).  However,  the  use  of  the 
envelope  equation  implies  that  the  wavelength  of  the  carrier  waves  is  much  smaller 
that  the  envelope  width,  which  often  seems  not  to  be  the  Ceise  in  visual  observations. 
Moreover,  as  pointed  out  by  Wu  (in  discussion  of  Akylas  et  al.  1988),  the  carrier 
waves  are  not  perpendicular  to  the  envelope  track  and  hence  it  is  an  open  question 
whether  the  asymptotic  radiation  of  the  envelope  wave  packets  can  be  satisfactorily 
explained  by  local  group  velocity  considerations. 

The  observed  wake  of  the  Coast  Guard  cutter  was  distorted  by  ambient  waves, 
thereby  limiting  the  comparison  to  averaged  properties  (Brown  et  al.  1989).  There 
was  still  variability  in  the  reduced  data  even  when  the  components  due  to  the  ambient 
waves  were  filtered.  This  could  be  due  to,  in  addition  to  the  ambient  waves,  variation 
in  ship  speed  or  heading;  time-dependent  ship  motion;  intersxtion  with  the  transverse 
Kelvin  wave  and  the  time-dependency  or  instability  of  the  solitary  wave  envelope. 
Other  sources  of  unsteadiness  are  the  propeller-induced  flow  field,  vortex  shedding 
(von  Karmen  vortex  street  type)  and  the  turbtilent  flow  in  the  neair  field.  Unsteady 
features  of  the  wake  may  play  a  very  important  role  in  the  generation  of  the  inner- 
angle  wavepackets.  It  is  therefore  desirable  to  investigate  the  unsteady  phenomena 
of  the  inner-angle  wavepackets. 
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Unsteady  waves  have  long  been  studied  in  seakeeping.  However,  the  main  objec¬ 
tive  in  seaJceeping  is  to  determine  the  hydrodynaunic  forces  acting  on  the  ship  and 
the  ship  motion  while  little  attention  is  given  to  the  fa^-field  wave  patterns.  Eggers 
(1957)  studied  the  wave  pattern  due  to  a  pulsating  translating  source  and  showed 
that  there  were  three  systems  of  constant  phase  curves.  Recent  work  by  Noblesse 
&  Hendrix  (1990)  showed  similar  results.  These  studies  have  shown  that  the  wave 
boundaries  change  with  the  frequen<y  of  the  oscillating  source.  Newman  (1961)  stud¬ 
ied  the  wave  pattern  due  to  a  non-oscillating  source  translating  in  an  incident  wave 
amd  showed  two  systems  of  “distorted  Kelvin  waves”.  All  these  studies  used  the 
method  of  stationary  phase.  Jankowski  (1990)  studied  linear  time-harmonic  wave 
patterns  caused  by  translating  and  pulsating  sources  using  fundamental  solutions 
satisfying  a  linear  free  surface  boundary  condition.  Nakos  and  Sclavounos  (1990) 
calculated  the  time-harmonic  wave  patterns  by  a  modified  Wigley  hull  model  trans¬ 
lating  and  oscillating  in  heave  using  a  R^kine  panel  method.  However,  none  of 
the  described  work  on  unsteady  wave  patterns  specifically  tried  to  explain  the  inner- 
«uigle  wavepacket  phenomenon.  Only  recently,  Mei  (1991)  used  a  unsteady,  nonlinear 
Schrodinger  equation  to  study  the  inner-angle  wave  phenomenon. 

Since  we  are  primarily  interested  in  far-field  waves,  we  simply  model  the  ship  by 
a  source-sink  pair  moving  below  the  free  surfawie  or  a  moving  pressure  distribution  on 
the  free  surface.  Unsteady  wakes  can  be  classified  into  two  types:  “diffraction”  waves 
generated  by  a  non-oscillating  disturbance  moving  in  incident  waves;  and  “radiation” 
waves  generated  by  an  oscillating  disturbance  moving  in  an  initially  undisturbed 
water.  Both  types  exist  for  am  advancing  ship  in  a  seaway,  since  the  seaway  is  not 
only  diffracted  by  the  presence  of  the  ship  but  also  causes  the  ship  motion,  which  in 
turn  generates  radiation  waves.  We  only  study  the  “radiation”  wake  here,  leaving 
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the  “diffraction”  type  and  the  interaction  between  the  two  types  of  wakes  for  future 
study. 

The  following  three  methods  are  used  in  our  study: 

1.  Nonlinear  wave  calculations:  It  is  rather  straightforward  to  apply  the  method 
described  and  used  in  the  previous  chapters.  However,  a  large  computational 
domain  on  the  free  surface  is  required  because  of  our  interests  in  the  far  field, 
resulting  in  a  large  number  of  unknowns  in  numerical  computations  (typically 
the  number  of  nodes  N  >  3000).  Some  special  techniques,  such  as  a  fast 
iterative  solver,  preconditioning  and  the  multigrid  method,  are  examined  «ind 
employed  to  speed  up  the  computations. 

2.  Linear  wave  calculations:  The  time  domain  Green  function  satisfying  the  linear 
free  surface  conditions  and  radiation  conditions  (Liapis  &  Beck  1986,  King 
1987,  Beck  &  Magee  1990)  is  used  to  calculate  the  wave  elevation. 

3.  Method  of  stationary  phase:  Although  it  is  not  easy  to  evaluate  wave  elevation 
by  the  method  of  stationary  phase,  the  wave  patterns  determined  by  constant 
phase  can  be  obtained  without  a  great  deal  of  computational  effort.  The  con¬ 
stant  phase  lines  are  obtained  in  a  similar  manner  as  Yih  &  Zhu  (1989a,b)  with 
some  modification.  In  Yih  &  Zhu  (1989a,b),  patterns  of  steady  ship  waves  are 
analyzed  with  the  added  effects  of  finite  depth,  surface  tension,  and  wake  ve¬ 
locity  profiles. 

6.1  Nonlinear  Wave  Calculations 

In  the  nonlinear  calculations  with  a  large  number  of  nodes  {N  >  3000),  am  itera¬ 
tive  solver  is  used  to  solve  the  system  of  linear  equations  representing  the  discretized 
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boundary  integral  equation: 

Ax  =  b.  (6.1) 

It  is  desirable  to  use  a  technique  to  accelerate  the  convergence.  Two  techniques,  the 
multigrid  method  and  matrix  preconditioning,  are  tested  with  two  iterative  bolvers, 
GMRES  and  a  successive  over- relaxation  algorithm  (SOR). 

6.1.1  Multigrid  and  matrix  preconditioning 

Multigrid  methods  have  been  widely  used  and  proved  effective  in  various  finite 
difference  schemes  (the  working  equation  being  of  lifferential  type).  Multigrid  meth¬ 
ods  have  been  infrequently  applied  to  boundary  integral  methods  (integral  type).  A 
brief  description  of  the  multigrid  method  used  with  the  GMRES  and  SOR  methods 
for  the  algebraic  system  resulting  from  the  desingularized  algorithm,  ar.d  the  dis¬ 
cussions  of  the  test  results  are  given  in  Appendix  A.  Unfortunately,  the  test  results 
show  that  the  multigrid  method  does  not  help  in  our  nonlinear  wave  calculations. 
The  GMRES  algorithm  without  a  multigrid  correction  is  superior. 

The  basic  idea  of  preconditioning  is  to  find  an  approximate  matrix  A  to  A  so 
that  A~^A  closely  approximates  the  identity  matrix,  yielding  a  strongly  diagonally- 
dominate  matrix.  The  preconditioned  equation; 

A-*Ax  =  A-*b  (6.2) 

is  then  better  conditioned  than  the  original  system  (6.1)  and  can  be  solved  with  much 
fewer  iterations.  In  general,  it  is  not  easy  to  find  a  good  approximate  A  matrix  for 
a  dense,  nonsymmetric  A  matrix.  However,  we  choose  A  as  A©,  t.e.  A  at  t  =  0. 
Although  the  free  surface  shape  ch2Uiges  with  time  and  so  does  A,  A,,  is  still  a  very 
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good  approximation  to  A,  especially  when  the  boundary  integral  is  desingularized. 
More  detail  on  the  test  results  for  the  nonlinear  wave  calculations  with  precondi¬ 
tioning  are  given  in  Appendix  B.  We  found  that  this  preconditioning  matrix  is  very 
effective  in  accelerating  the  convergence.  For  example,  preconditioning  saves  over 
80%  of  the  total  CPU  time  to  solve  the  system  of  3367  linear  equations  in  the  wave 
computation  for  1000  time  steps  shown  in  Fig.  6.1.  The  preconditioning  technique, 
not  the  multigrid  method,  is  therefore  used  in  the  nonlinear  wave  computations. 

6.1.2  Computational  results 

We  first  calculate  the  waves  generated  by  a  submerged  source-sink  pair.  The 
separation,  the  depth  of  the  submergence,  and  the  strength  of  the  source  and  the 
sink  are  chosen  to  approximately  model  the  Coeist  Guard  cutter.  The  cutter  h2is  a 
25.3m  waterline  length,  5.2m  beam,  1.8m  draught,  and  the  ship  speed  U  was  7.7m/ s 
(Brown  et  al.  1989).  The  source-sink  pair  starts  to  move  horizontally  from  the  rest. 
We  use  the  same  separation  (87.7%  of  the  waterline  length),  submerged  depth  (50% 
of  the  draught)  and  strength  {UA^,  where  A„  =  3.3m*  is  the  submerged  cross- 
sectional  area  near  midship)  of  the  source-sink  pair  as  used  by  Hall  &  Buchsbaum 
(1990).  The  problem  is  made  dimensionless  based  on  the  ship’s  waterline  length.  The 
strength  of  the  pair  are  multiplied  by  (1  -  e~'‘‘)  where  fi=2.0.  The  time  simulation 
is  sufficiently  long  so  that  the  waves  behind  the  pair  achieve  a  steady  state.  Fig.  6.1 
shows  contour  and  isometric  plots  of  the  waves  when  the  pair  has  traveled  about  16 
ship  lengths.  The  wake  is  essentially  a  Kelvin  wake  although  there  is  interference 
between  the  wave  systems  by  the  source  and  the  sink.  This  interference  does  not 
form  a  wavepacket  like  the  one  observed  in  the  experiment  of  the  cutter  Point  Brower 
at  an  inner-angle  of  less  than  19.5  degrees. 
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Fig.  6.2  shows  the  waves  generated  by  two  moving  dipoles.  The  dipole  strength 
is  made  as  large  as  possible  so  that  the  waves  are  nonlinear  but  not  breaking.  The 
submergence  depth  of  the  dipole  is  set  to  1,  the  separation  of  the  dipoles  is  10,  and 
the  strength  is  2.5.  The  depth  Froude  number  is  1.0.  Fig.  6.3  and  Fig.  6.4  show  the 
waves  due  to  two  moving  free  surface  pressure  patches.  Each  patch  has  a  length  1 
and  beam  0.2.  The  pressure  distribution  for  each  patch  is 


where 


1  |x|  <  0.8 

/»*(x)  =  <  i{cos  [57r(i  ^  0.8)]  +  1}  0.8  <  ±x  <  1 

0  lz|  >  1.0 


1  |y|  <  0.1 

=  ‘  i{cos[107r(y:F  0.1)1 +  1}  0.1  <  ±y  <  0.2 

0  |y|  >  0.2 

V 

and  (i,  y)  are  the  local  coordinates  attached  to  the  centroid  of  a  patch.  The  centroids 
of  the  two  pressure  patches  are  separated  by  1.8.  In  Fig.  6.3,  the  pressure  is  small 
{po  =  0.001)  so  the  waves  are  essentially  linear.  In  Fig.  6.4,  the  pressure  is  20  times 
larger  so  the  waves  are  nonlinear  (the  maximum  wave  slope  is  about  0.30).  Fig.  6.5 
compares  the  longitudinal  wave  cuts  of  the  linear  and  nonlinear  steady  waves  at 
different  y  when  the  linear  waves  of  Fig.  6.3  are  multiplied  by  20.  As  seen  in  Fig.  6.5, 
the  nonlinear  effect  is  not  very  significant  and  again  no  small  angle  wavepackets  form. 
Instead,  the  two  wave  systems  by  the  forv  ird  and  aft  disturbances  seem  parallel  to 
each  other.  The  computations  (not  shown  here)  using  other  separation  distances  of 
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the  disturbances  reveal  the  same  feature,  except  when  the  separation  distance  is  so 
small  that  the  two  wave  systems  merge  into  one  Kelvin-type  wave  system. 

We  then  calculate  the  unsteady  wakes  generated  by  one  pressure  patch  of  form 
(6.3)  translating  at  a  speed  and  pulsating  with  a  frequency  u  such  that  po  = 
Po  sinw.  Figures  6.6  to  6.11  show  the  waves  caused  at  diflferent  vadues  of  the  reduced 
frequency  r  =  with  =  0.4.  For  small  t,  two  wave  systems  can  be  seen, 

one  with  a  bounding  angle  Izu-ger  than  the  Kelvin  angle  and  the  other  smaller  than 
the  Kelvin  angle.  As  the  frequency  increases,  only  one  wave  system  is  visible  with 
a  smaller  bounding  angle.  Although  it  is  not  easy  to  measure  bounding  angles 
precisely  from  these  contour  plots,  it  is  cleau^  that  the  angles  decrease  with  increasing 
disturbance  frequency. 

Since  the  pressure  patch  is  not  a  point  pressure,  a  Proude  number  effect  is  ex¬ 
pected.  Fig.  6.12  shows  wave  contour  plots  for  =  0.2  and  or  =  5  keeping  r  =  1.0. 
The  wave  orientation  in  relation  to  the  wavepacket  is  quite  different  compared  to 
that  in  Fig.  6.8  (with  the  same  vadue  of  t  but  Fn  =  0.4)  although  the  bounding 
packet  angle  is  not  affected. 

6.2  Linear  Wave  Calculations 

The  steady  nonlinear  wave  computations  have  shown  that  nonlinear  effects  are 
not  essentiad  to  the  generation  of  the  inner-amgle  wavepackets  while  the  unsteadiness 
of  the  flow  can  cause  the  inner-angle  waves.  Further  study  on  the  features  of  the  un¬ 
steady  wake  using  linear  computations  is  suggested.  We  will  now  study  the  unsteady 
wakes  by  linear  computations  using  a  time-domain  Green  function. 
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6.2.1  Time-domain  Green  function 

The  time-domain  Green  function  G{P^P\t,V)  we  use  is  the  potential  due  to  an 
impulsive  source  at  time  t'  and  location  P’: 

G  (P,  P', t')  =  (;  -  5 (<  -  0  +  G (P, P\t,l!)H{t- 1>) .  (6.4) 

Here, 

G  (P,  P\  t,  i')  =  j^dkyJTg  sin  (t  -  t'))  e*<*+*'>  Jo  (kR)  (6.5) 

is  the  wavy  part,  P  =  (i,  y,z)  is  the  field  point,  and  P'  =  (i'(f^),y'(<'),z'(f'))  is  the 
source  point,  r  =  [(x  —  x')^  +  (y  ~  ~  ^0^1*  's  the  distance  between  the  field 

and  source  points,  f  =  [(x  —  x')^  +  (y  “  yO^  +  (2  +  is  the  distance  between  the 
field  point  and  the  image  of  the  source  point,  and  P  =  [(x  —  x')*  -f  (y  —  J..  is 

the  zero-order  Bessel  function,  S{i)  is  the  delta  function,  and  H{t)  is  the  unit  step 
function. 


The  Green  function  G{P,P*,t,t*)  satisfies  the  following  governing  equation, 

V*G  =  -Air6iP-F)6(t-t'), 

with  the  linear  free  surface  condition  on  x  =  0, 

d^G  dG 

an+i'-fc  =  “  (onr  =  0), 

and  the  initi2J  and  far-field  conditions, 

^  dG  ^ 

G,  -^  =  0  t  -  t'  <  0 
VG  —>0  as  r  — f  00 


The  potential  due  to  a  source  of  strength  a{t)  moving  in  the  path,  (x'  =  x'(f),y'  = 
y'{t),z'  =  z\i))  is  then 

^(P,f)=  f*  a{t')GiP,P\t,f)di\ 

JO 


(6.6) 
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The  free  surface  elevation  in  the  fixed  coordinate  system  is, 

1 

V{x,y,t)=^ — (on  z  =  0).  (6.7) 

g  at 

When  X  is  replaced  hy  x  —  Vt  and  ^  by  ^  one  obtains  the  wave  elevation  in 

the  coordinate  system  moving  with  velocity  V  in  the  —  x  direction, 

J7(i,y,<)  = -j(^  + (on  r  =  0).  (6.8) 

We  use  the  subroutine  initially  developed  by  Liapis  (1986)  and  later  improved 
by  King  (1987)  and  Magee  (1990)  to  calculate  G(P,  P',<,t').  The  time  convolution 
integral  in  (6.6)  is  performed  using  a  simple  trapezoidal  nile. 

6.2.2  Computational  results 

Since  discretization  of  the  free  surface  is  not  necessary  in  the  linear  calcula¬ 
tions,  less  computation  time  is  required  compared  to  the  nonlinear  wave  cadculations. 
Therefore,  we  are  able  to  consider  a  larger  free  surface  area  with  a  finer  grid.  The 
disturbance  used  is  a  source-sink  pair.  The  submergence  depth  is  1.  The  separation 
of  the  source  and  sink  is  also  1.  The  problem  is  made  dimensionless  based  on  the 
depth.  Figures  6.13  to  6.17  show  the  wave  patterns  by  a  translating  and  pulsating 
source-sink  pair  with  zero-mean  strength  for  different  vadues  of  t.  Similar  to  the  non¬ 
linear  computations,  two  wave  systems  can  be  identified  at  low  frequency.  Again, 
as  the  frequency  increases,  the  angles  of  the  two  wave  systems  merge  and  only  one 
system  is  seen.  The  angles  measxrred  from  these  figures  are  shown  in  Fig.  6.18  as  a 
function  of  r.  The  inner  wave  angle  approaches  the  Kelvin  amgle  as  r  approaches 
zero.  The  outer  wave  angle  has  a  maximum  around  90  degrees  at  r  =  1/4  and  de¬ 
creases  to  the  Kelvin  angle  as  t  approaches  zero.  From  Fig.  6.18,  the  inner-angle  of 
the  Point  Brower  observations  corresponds  to  1.0  <  t  <  1.5. 
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Fig.  6.19  shows  contour  and  isometric  plots  of  the  far-field  waves  generated  by 
a  translating  and  pulsating  source-sink  pair  with  nonzero  mean  strength  which  may 
correspond  to  the  surge  motion  of  a  ship.  The  puslation  amplitude  is  the  same  as 
the  mean  value.  The  reduced  frequency  is  r  =  1.0.  An  inner-angle  wavepacket  at 
about  14  degrees  is  observed. 

Figures  6.20  (t  =  1.0)  and  6.21  (t  =  2.0)  show  the  wave  patterns  corresponding 
to  those  in  Figs.  6.16  and  6.17,  respectively,  with  the  Froude  number  doubled 
and  the  frequency  u  halved  to  keep  r  unchanged.  The  computational  domain  on  the 
free  surface  is  enlarged  by  as  explained  by  (6.25)  in  the  following  section.  The 
wave  systems  with  larger  bounding  angles  in  Fig.  6.16  and  Fig.  6.17  are  not  visible 
in  Fig.  6.20  and  Fig.  6.21  (the  wave  elevation  may  be  too  small  to  be  caught  in  the 
contour  plots).  Again,  similar  to  the  nonlinear  calculations  (Fig.  6.8  and  Fig.  6.12), 
the  wave  orientations  along  the  wavepackets  in  Fig.  6.20  and  Fig.  6.21  are  slightly 
different  from  those  in  Fig.  6.16  and  Fig.  6.17  while  the  packet  angles  are  almost 
unaffected. 

6.3  Wave  Patterns  by  Method  of  Stationary  Phase 

The  method  of  stationary  phase  provides  an  easy  way  to  determine  far-held  wave 
patterns  in  terms  of  curves  of  constant  phase.  The  disturbance  can  be  a  ship,  a  free 
surface  pressure  distribution,  underwater  body  or  singularity  (source-sink  pair).  In 
this  analysis,  the  particularity  of  the  disturbance  does  not  affect  the  establishment  of 
the  formulas  for  determining  the  phase  lines  (Yih  &  Zhu  1989a, b).  If  the  region  under 
consideration  is  sufficiently  far  from  the  disturbance,  the  wave  pattern  obtained  may 
be  independent  of  the  local  disturbance  details. 
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6.3.1  Constant  phase  lines 

The  velocity  potential  $(x,  t)  can  be  expressed  as 


«(£,()  =  iee{«(f)e“}, 


(6.9) 


assuming  the  system  is  linear  and  has  reached  a  harmonic  limit  (^cle.  Clearly, 
the  spatial  function  <i>  also  satisfies  Laplace  equation.  The  linearized  free  surface 
boundary  condition  for  $  in  the  coordinate  system  moving  with  the  disturbance  is 


Substituting  (6.9)  into  (6.10)  gives  the  free  surface  condition  for  <f>, 

+  =  0  (ODZ  =  0). 

Separation  of  variables  gives 


<i>  oc  c**cxp  [t  {kgX  +  fcyy)] 


(6.10) 


(6.11) 


(6.12) 


where  the  wavenumber  vector  is  given  by 


{kx,ky)  =  A:(cos0,sin^). 


(6.13) 


Substitution  of  (6.13)  and  (6.12)  into  (6.11)  gives  the  relation  between  k  and  6, 


k  =  {kFn  cos  9  +  u;)*. 

Solving  (6.14)  for  k  gives 

.  _  1  +  2t  cos  9  ±  \/l  +  4r  cos  9 
2Fn^  cos^  9 

with  its  first  derivative 

^^-..-.(Vl+4rcos»±l) 
d9  y/l  -f-  4t  cos  9 


(6.14) 


(6.15) 


1 


(6.16) 
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where  r  =  F„u;  is  the  reduced  frequency.  Since  k  is  real,  1  +  4t  cos  ^  >  0,  the  region 
for  &  must  be  restricted  to  0  <  e<e,,  where 

f  IT  if  T  <  0.25 


^0  = 


[  JT  —  COS~*  (^)  if  T>  0.25 


The  spatial  potential  <f>  can  then  be  expressed  as 

^(x)  =  f  d6A{6)e’‘^exp\i  {kxX kyy)],  (6.17) 

Jo 

where  A{0)  is  a  function  depending  on  the  particular  disturbance.  The  wave  elevation 
f}(x,t)  has  a  form  similar  to 

f}{x,t)  =  i2e{»7(x)c“*'‘},  (6.18) 

where  J7(x)  has  a  form  similax  <{>  (6.17), 

ri{x)=  d^B(5)c‘ <**•+"*»>  =  (6.19) 

Jo  Jo 

and  B{d)  is  a  function  depending  on  the  the  particular  disturbance,  R  and  /?  are  the 
polar  coordinates  of  the  field  point. 


(i,  y)  =  R(cos  0,  sin  ^) 


and 


rf)  =  k(cos  0  cos  5  +  sin  sin  0). 


Evaluation  of  (6.17)  by  the  stationary  phase  method  (Copson  1965)  gives 


M 


1/2 


^  (g  20) 


[2Rr(M)\ 

where  the  summation  is  made  over  all  the  stationary  points  5,-  which  cire  given  by 
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or 

X  dkg  +  y  dky  =  0.  (6.21) 

The  “±”  sign  in  (6.20)  is  the  same  as  the  sign  of  the  second  derivative  V>"(/S,  9i).  For 
simplicity,  the  subscript  t  and  the  baur  over  6  are  omitted  hereafter. 

The  curve  of  constant  phase  can  be  obtained  without  the  knowledge  of  A{9)  and 
B{9).  The  curve  of  constant  phase  C  is  given  by 
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The  sign  of  c  must  be  properly  chosen  to  satisfy  the  radiation  conditions  (lines 
of  constant  phcise  do  not  extend  to  x  — +  — oo).  At  the  first  glance,  there  seems  to 
be  two  sets  of  constant  phase  lines  since  (6.14)  hats  two  roots  for  k.  By  carefully 
examining  (6.25),  we  find  that  (6.25)  using  positive  c  and  the  sign  satisfies  the 
radiation  condition  for  0®  <  5  <  90"  but  not  for  90"  <  $  <  9o.  The  phase  c  must 
change  sign  for  90"  <  6  <  So  (6.25)  using  the  set  gives  two  systems  of  phase 
lines.  For  the  ”  sign  in  (6.25),  c  must  be  positive  so  that  the  radiation  condition 
is  satisfied.  Equation  (6.25)  using  ”  gives  only  one  system  of  constant  phase  lines. 

Figures  6.22  to  6.29  show  lines  of  constant  phase  C  for  several  values  of  t.  Each 
figure,  for  a  given  t  >  0,  contains  three  distinct  systems  of  waves  similar  to  those 
in  Noblesse  &  Hendrix  (1990)  using  a  different  approach.  For  r  =  0  (steady  waves), 
the  ”  set  of  the  constant  phase  lines  disappear  (the  wavelength  becomes  infinitely 
large),  while  the  first  and  second  systems  of  the  phase  lines  become  identical  to  the 
normal  Kelvin  wave  pattern.  The  first  system  is  very  similar  to  the  Kelvin  wave 
pattern  which  consists  of  transverse  waves  and  diverging  waves  auid  has  a  cusp  line 
whose  angle  is  less  than  the  Kelvin  angle  19.5".  The  second  system  also  consists 
of  transverse  waves  and  diverging  waves  for  t  <  0.25  with  a  cusp  line  of  a  larger 
angle  than  19.5".  For  t  >  0.25  the  transverse  waves  disappear.  The  third  system 
has  ring-type  waves  for  t  <  0.25.  For  t  >  0.25,  part  of  the  ring-type  waves  become 
diverging  waves. 

A  line  of  constant  phase  C  is  disconnected  at  the  cusp  line  of  a  wave  system 
since  the  constant  c  has  a  jump  of  7r/2  across  the  point.  This  is  because  the  second 
derivative  of  changes  its  sign  across  the  point,  which  will  be  explained  in  more 

detail  in  the  following  sections.  Because  of  this  jump,  the  crests  (or  trough)  of  the 
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transverse  waves  and  the  crests  (or  trough)  of  the  diverging  waves  do  not  meet  on 
the  cusp  line.  The  phase  lines  shown  in  Noblesse  &  Hendrix  (1990)  correspond  to 
the  lines  of  constant  c  in  our  case  where  the  lines  should  be  continuous. 


6.3.2  Bounding  angles  for  constant  phase  lines 


From  (6.25),  we  have,  for  the  first  and  second  systems. 


|y| 


sin  9  cos  0 

sin^  9  +  v^TTTrcosd’ 


(6.26) 


where 


/l,2  = 


1  -f  2t  cos  9  +  y/1  +  cos  9 
1  +  2r  cos  9  +  y/1  +  4t  cos  9 


and  for  the  third  system. 


|y| 


sin  d  cos  d 

sin*  0  —  >/l  +  4t  cos  0  ’ 


(6.27) 


where 


/a  = 


1  +  2r  cos  9  —  y/1  +  At  cos  0 
1  +  2r  cos  9  —  \/r'+4Tcosl0 


Equation  (6.26)  is  plotted  in  Fig.  6.30  for  various  values  of  t.  The  curves  for  0®  < 
9  <  90®  correspond  to  the  first  syst«n  of  phase  lines  and  9  >  90®  corresponds  to  the 
second  system  of  phase  lines  as  explained  above.  The  third  system  (6.27)  is  plotted 


in  Fig.  6.31  for  various  values  of  r. 


The  bounding  2uigle  for  each  system  can  be  determined  from  Fig.  6.30  and 
Fig.  6.31  or  from  (6.26)  and  (6.27)  or  (6.25)  directly.  The  maximum  of  \y\/x  of  each 
system  gives  the  bounding  angle  a  =  tan-^(|y|/x).  For  the  first  system,  0®  <  ^  <  90® 
in  Fig.  6.30,  |y|/x  has  one  maximum  which  is  also  a  local  maximum  (*.e.  the  deriva¬ 
tive  of  |y|/x  is  zero).  For  the  second  system,  90®  <  <  9^  in  Fig.  6.30,  ly|/i  also 
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has  one  maximum  which,  however,  is  not  a  local  maximum  when  r  >  1/4.  For  the 
third  system.  Fig.  6.31,  ly|/x  can  have  a  maximum  of  ±oo  for  r  less  than  a  certain 
value.  |y|/x  heis  a  jump  from  oo  to  — oo  (or  — oo  to  oo)  across  an  angle  Be  which 
is  the  root  of  sin*0  —  y/\  +  4t  cos  9  =  0.  The  infinity  of  |y|/x  does  not  correspond 
to  the  angle  of  the  bounding  line.  The  angle  of  the  bounding  line  in  this  system 
only  comes  from  the  local  maximum  of  |y|/x  where  the  derivative  of  |yl/x  is  zero. 
Fig.  6.32  shows  tan“^(|y  |/x)  as  a  function  of  6  for  the  third  system.  From  Fig.  6.32,  it 
is  easier  to  determine  the  angle  of  the  boimding  line  by  finding  the  global  maximum 
of  tan“^(|y|/x). 

Fig.  6.33  shows  the  bounding  angles  of  the  three  systems  (denoted  by  oj,  and 
03  respectively)  as  a  function  of  t.  oi  is  a  continuous  function  of  t.  oj  has  a  jump 
from  54.74"  to  90"  and  03  has  a  jump  from  180"  to  126.26"  when  t  changes  from 
0.25—  to  0.25+.  Fig.  6.33  is  convenient  to  determine  the  angles  of  the  inner-angle 
wavepackets.  Fig.  6.30  and  Fig.  6.31  help  estimate  the  decay  rates  of  the  waves. 

6.3.3  Decay  rates  of  waves 

The  wave  decay  rate  for  large  R  can  be  estimated  by  applying  the  method  of 
stationary  phase  to  (6.19)  without  knowing  B{9)  (Copson  1965).  The  decay  rate  is 
closely  related  to  the  second  derivative  of  x(?  (corresponding  to  the  first  derivative 
of  |y|/x,  Newman  1977).  Fig.  6.30  and  Fig.  6.31  are  helpful  since  they  provide  the 
information  about  the  derivative  of  |y|/x. 

•  Along  /S,  in  which  ^  ^  o,-  and  P  <  max(o,),  i  =  1,2,3,  the  wave  elevation 
is  given  by  (6.20)  since  the  second  derivative  of  ^  (corresponding  to  the  first 
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derivative  of  |y|/a;,  Newman  1977)  does  not  vanish  at  the  stationary  points. 
The  waves  decay  at  a  rate  of  For  simplicity,  we  express  this  decay  zis 

,(/?,« oc^  +  0(i).  (6.28) 

•  If  is  greater  than  max(ai),  there  are  no  stationary  points  and  the  waves  decay 
faster, 

oc  ;!  +  (6-29) 

•  The  decay  rates  along  the  bounding  angles  of  the  wave  systems  axe  as  follows: 

—  Along  =  Oj,  for  all  r,  is  zero  for  the  first  wave  system  at  the  sta¬ 
tionary  point  while  0"  is  not  zero  for  the  other  two  wave  systems.  The 
contribution  from  the  first  system  dominates  and  decays  as  while 

the  contributions  from  the  other  two  systems  decay  as  Therefore, 

the  decay  rate  along  =  oj  is 

r}{R,0)  oc  d"  ^(■^1/2  )•  (6.30) 

—  Along  0  —  02,  is  zero  for  the  second  system  at  the  stationary  point 
only  when  r  <  0.25.  When  t  >  0.25,  is  not  zero.  Also,  is  also  not 
zero  for  the  other  two  systems  for  all  t.  Therefore,  the  decay  rate  along 
02  is 

if  T  <  0.25 

ri{R,0)  oc  < 

;R^d-0(^)  if  T>  0.25. 

-  Along  0  =  03,  is  zero  for  the  third  system  for  r  >  0.25.  0"  is  not  zero 
for  the  other  two  wave  systems  for  all  r.  Therefore,  the  waves  decay  as 
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^  +  0(i)  if  r  >  0.25 

77(/2,/9)  a  - 

+  ^(«)  if  T  <  0.25. 

6.4  Comparisons  and  Discussions 

We  now  compare  the  results  of  the  nonlinear  ’^nd  linear  calculations  as  well  as 
the  stationary  phase  analysis,  and  steady  waves  to  unsteady  waves. 

Fig.  6.34  shows  the  contour  plot  of  the  waves  generated  by  the  same  source-sink 
pair  used  in  Fig.  6.1  to  simulate  the  cutter  Point  Brower  using  the  linear  calculation. 
Both  the  nonlinear  (Fig.  6.1)  and  linear  (Fig.  6.34)  calculations  show  saime  interfer¬ 
ence  between  the  wave  system  caused  by  the  source  and  that  caused  by  the  sink. 
The  interference  is  not  clearly  separated  from  the  Kelvin  diverging  waves.  Fig.  6.35 
and  Fig.  6.36  show  results  of  the  linear  and  nonlinear  wave  calculations  for  the  same 
source-sink  pair  used  in  Fig.  6.1  and  Fig.  6.34.  Here  the  pair  not  only  translates  but 
also  heaves  about  the  mean  submergence  depth  used  for  Fig.  6.1.  We  use  a  heave 
amplitude  of  25%  of  the  mean  depth  and  a  reduced  frequency  t  =  0.5  although  the 
information  on  the  motion  of  the  cutter  is  not  known.  An  inner- angle  wavepacket 
is  seen  separated  from  the  Kelvin  diverging  waves  at  an  amgle  between  13.5®  and 
14.5®  in  both  linear  and  nonlinear  calculations.  The  wavepacket  by  the  nonlinear 
calculation  seems  more  isolated. 

Three  wave  systems  are  identified  in  the  stationary  phase  analysis.  However,  the 
phase  lines  do  not  give  the  quantitative  information  about  the  wave  elevation.  To 
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calculate  the  wave  elevation,  it  is  necessary  to  know  A{9),  or  which  depends  on 

the  psurticular  disturbance.  It  is  possible  that  one  may  not  be  able  to  identify  all  the 
three  wave  systems  in  the  isometric  and  contour  plots  of  the  wave  elevation.  At  low 
frequency,  we  see  the  first  and  second  systems  in  the  plots  of  the  linear  and  nonlinear 
waves.  The  third  system  is  not  seen  most  likely  because  it  has  a  large  wavelength, 
a  small  amplitude,  or  is  outside  the  computational  window.  At  an  intermediate 
frequency,  the  third  system  can  be  identified  while  it  is  difficult  to  distinguish  the 
second  system  from  the  first  system  because  the  bounding  lines  of  the  two  systems 
are  close  and  the  waves  interact.  At  high  frequency,  the  second  system  has  very 
short  wavelengths,  finer  than  our  computational  grids  can  resolve.  We  can  see  the 
third  system  or  the  combination  of  the  first  and  third  systems  but  it  is  difi&cult  to 
distinguish  them  from  the  plots  of  wave  elevation.  Fig.  6.37  shows  a  contour  plot 
of  wave  elevation  due  to  a  translating  pulsating  source-sink  pair  with  a  zero-mean 
strength  and  r  =  1.0  using  the  linear  calculation.  The  corresponding  constant  phase 
lines  are  overlayed  on  the  contour  lines.  The  agreements  in  the  wavelength,  wave 
orientations,  and  the  phase  in  the  inner  region  eind  outer  region  are  fairly  good. 

Figiires  6.38  to  6.42  show  the  wave  elevation  (shown  by  darkness  of  color)  due 
to  the  same  source-sink  pair  used  in  Fig.  6.37  but  with  t  =  0.5  at  five  different  times 
within  em  oscillation  cycle,  beginning  at  time  t  =  12.0  when  the  pair  has  traveled  18 
ship  lengths.  The  period  of  the  cycle  is  100  A  f  where  At  =  0.02  is  the  time  step  used 
in  the  calculation.  As  can  be  seen,  the  initial  transient  effects  have  died  out  and  the 
waves  have  almost  achieved  the  limit  cycle  since  Fig.  6.38  and  Fig.  6.42  2ure  nearly 
identical.  To  examine  how  these  unsteady  waves  move,  the  wave  elevation  within 
one  cycle  along  the  symmetry  plane  and  the  ray  of  the  inner-angle  wavepacket  of 
14  degrees  is  shown  in  Fig.  43  and  Fig.  44  respectively  (the  wave  motions  ran  be 
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seen  more  clearly  from  the  animation  of  the  waves).  In  Fig.  43  and  Fig.  44,  R  is 
the  distance  from  the  center  of  the  source-sink  pair.  As  expected,  the  transverse 
waves  (Fig.  6.43)  radiate  from  the  disturbance.  The  waves  along  the  line  of  14“ 
(Fig.  6.44)  also  radiate  from  the  disturbance  but  with  a  fluctuation  at  lower  frequency. 
The  animation  clearly  shows  that  the  inner-angle  wavepacket  originates  from  the 
disturbance  and  moves  away  from  it.  The  wavepacket  fluctuates  along  the  14“  line, 
which  agrees  with  the  observed  sinuous  fluctuation  along  the  “solitary”  feature  in 
the  wake  of  the  cutter  Point  Brower  (Brown  et  al.  1989).  The  animation  also  shows 
that  there  are  long  crest  waves  moving  away  from  the  wavepjuiket  periodically.  As  a 
wave  {e.g.,  a  wave  crest)  in  the  wavepacket  moves  along  the  14“  line,  half  of  it  merges 
into  a  long  crest  wave  and  the  other  half  keeps  moving  along  the  wavepacket. 

From  Fig.  6.43  and  Fig.  6.44,  it  can  be  seen  that  the  wave  amplitude  along  the 
wavepacket  is  larger  and  decays  at  a  slower  rate  than  that  along  the  symmetry  plane. 
Moreover,  Fig.  6.45  shows  the  peak  wave  amplitude  within  the  diverging  wave  system 
of  the  steady  wake  by  the  nonlinear  calculation  shown  in  Fig.  6.1.  The  peak  wave 
amplitude  at  a  distance  x  is  obtained  by  seeching  the  maximum  wave  amplitude 
in  the  y  direction  within  the  diverging  wave  region  (15“  <  /3  <  19.5“).  Fig.  6.46 
shows  the  peak  wave  amplitude  within  the  region  (10®  <  l3  <  15“)  of  the  interference 
between  the  bow  and  stern  wave  systems  in  Fig.  6.1.  Fig.  6.47  2md  Fig.  6.48  show 
the  peak  wave  amplitudes  within  the  diverging  wave  region  (15“  <  /?  <  19.5“)  and 
interference  region  (10“  <  <  15“)  for  the  steady  wake  by  the  linear  calculation 

shown  in  Fig.  6.34.  Fig.  6.49  shows  the  pe2d«  amplitudes  at  three  different  times 
within  the  the  inner-angle  wavepacket  region  (10“  <  /?  <  15“)  in  the  unsteady 
wake  by  the  nonlinear  calculation  shown  in  Fig.  6.36.  Fig.  6.50  shows  the  the  peak 
amplitudes  at  the  five  different  times  in  one  oscillation  cycle  within  the  the  inner- 
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angle  wavepacket  region  in  the  unsteady  waJce  by  the  linear  calculation  shown  in 
figures  6.38  to  6.42. 

The  upper  envelope  of  the  amplitude  in  each  figures  (figures  6.45  to  6.50)  gives 
the  decay  rate  of  the  wavepacket  amplitude.  A  straight  line  showing  the  algebraic 
decay  rate  by  the  stationary  phase  analysis  is  drawn  at  the  upper  bound  of  the  peak 
2unplitude  in  each  figures.  As  seen  from  these  figures,  for  the  steady  wake,  the  waves 
by  both  nonlinear  and  linear  calculations  decay  at  a  rate  close  to  within 

the  diverging  wave  region  while  at  within  the  interference  region  in  the 

steady  wake.  This  is  consistent  with  the  results  of  the  normal  steady  Kelvin  wake 
theory  (Ursell  1960).  For  the  unsteady  wake,  the  waves  by  both  nonlinear  and  linear 
calculations  decay  at  a  rate  close  to  0{R~^^^)  within  the  inner-angle  wavepacket, 
which  is  consistent  with  the  result  of  the  stationziry  phaise  analysis  in  this  thesis. 
The  experimental  data  of  Brown  et  al.  (1989)  are  also  redrawn  in  Fig.  6.51  to  show 
the  algebraic  decay  of  the  peak  amplitude.  The  observed  decay  rate  (Fig.  6.51)  is 
closer  to  0{R~^^^)  of  the  unsteady  inner-angle  wavepacket.  Brown  et  aL  (1989) 
compare  the  experimental  data  to  the  0{R~^^^)  decay  of  a  steady  linear  Kelvin  wake 
theory  to  draw  the  conclusion  that  nonlinear  effects  play  an  important  role  in  the 
generation  of  the  inner-angle  wavepacket. 

6.5  Conclusions 

Our  study  shows  that  nonlinear  effects  are  not  necessary  for  the  generation  of 
inner-angle  wavepackets;  unsteady  effects  can  cause  the  inner-angle  waves.  The 
observed  inner-angle  of  10.8  degrees  in  the  full  scaJe  measurement  on  the  cutter  Point 
Brower  may  correspond  to  a  ship  motion  at  a  reduced  frequency  within  1.0  <  t  <  1.5. 
In  the  unsteady  wake,  the  wave  amplitude  decays  as  along  the  ray  of  the 
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inner  packet,  the  same  rate  at  which  the  normal  steady  Kelvin  wave  decays  along 
its  cusp  line  (Ursell  1960).  Hence,  it  appears  that  the  persistence  of  the  inner-angle 
wavepacket  behind  the  ship  Ccin  be  explained  by  unsteady  linear  wave  theory  alone. 
Our  study  also  shows  that  waves  can  be  confined  in  a  very  narrow  region  if  the 
frequency  is  very  high.  Therefore,  the  high  frequency  propeller-induced  unsteady 
flow  may  also  contribute  to  the  very  narrow  F  structure  in  the  wake. 

The  Froude  number,  based  on  the  pressure  patch  length  in  our  nonlinear  calcu¬ 
lations  or  the  submergence  depth  of  the  source-sink  pair  in  our  linear  wave  calcu¬ 
lations,  affects  the  wave  orientation  in  relation  to  the  wavepacket  without  affecting 
the  wavepacket  angles  if  the  reduced  frequency  remsuns  unchanged. 


Figure  6.1:  Nonlinear  waves  by  a  translating,  non-oscillating  source-sink  pair, 
F„=0.49,  free  surface  domain:  (0<i<ll,  0<y<  4.5)  Contour 
lines  are  apart  by  0.005 


Figure  6  2-  Nonlinear  waves  by  two  translating,  non-oscillating  dipoles,  Fn— 1.0,  free 
surface  domain:(0  <  i  <  4.5  0  <  y  <  1.8  ),  Contour  Unes  are  apart 

by  0.001 


Figure  6.3:  Nonlinear  waves  by  two  translating,  non-oscillating  surface  pressure 
patches  (weak  disturbance),  Contour  lines  are  apeurt  by  0.0002 

(0  <  X  <  5.2,0  <y<  1.6) 


Figure  6.4:  Nonlinear  waves  by  tvro  translating,  non-oscillating  surface  pressure 
patches  (strong  disturbance),  Contour  lines  are  apart  by  0.004) 

(0<x<5.2,0<y^l.6) 


Wave  elevation 
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Figure  6.5:  Comparison  of  Nonlinear  waves  to  “linear"  waves  by  two  translating, 
non-oscillating  surface  pressure  patches  (strong  disturbance) 


Figure  6.6:  Nonlinear  waves  by  a  translating  and  oscillating  surface  pressure  patch 
with  zero  mean,  (t  =  0.25,  Fn  =  0.4),  Contour  lines  are  apart  by  0.0001 

(0  <x  <5.2,0  <y<  1.6) 


Figure  6.7:  Nonlinew  waves  by  a  translating  and  oscillating  surface  pressure  patch 
with  zero  mean,  (t  =  0.5,  Fn  =  0.4),  Contour  lines  are  apart  by  0.0001 

(0<a:<5.2,0<y  <  1.6) 
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Figure  6.8:  Nonlinear  waves  by  a  translating  and  oscillating  surface  pressure  patch 
with  zero  mean,  (t  =  1.0,  F’„  =  0.4),  Contour  lines  are  apart  by  0.0001 

(0<x<5.2,0<y<1.6) 


Figure  6.9:  Nonlinear  waves  by  a  translating  and  oscillating  surface  pressure  patch 
with  zero  mean,  (r  =  2.0,  F„  =  0.4),  Contour  lines  are  apart  by  0.00005 

(0<x<5.2,0<y<1.6) 


Figure  6.10:  Nonlinear  waves  by  a  translating  and  oscillating  surface  pressure  patch 
with  zero  mean,  (t  =  4.0,  =  0.4),  Contour  lines  are  apart  by  0.00005 

(0<x  <5.2,0  <y<  1.6) 


Figure  6.11:  Nonlinear  waves  by  a  translating  and  oscillating  surface  pressure  patch 
with  zero  mean,  (t  =  6.0,  Fn  =  0.4),  Contour  lines  are  apart  by  0.00001 

(0<x  <5.2,0  <y  <1.6) 


Figure  6.12:  Nonlinear  waves  by  a  translating  and  oscillating  surface  pressure  patch 
with  zero  mean,  (r  =  1.0,  F„  =  0.2),  Contour  lines  are  apart  by  0.0001 

(0  <x  <5.2,0  <y<  1.6) 
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Figure  6.13;  Linear  waves  by  a  translating  and  oscillating  source-sink  pair  with  zero- 
mean  strength,  (t  =  0.125)  (0<x<ll,0<y<  4.5) 


I 


Figure  6.14:  Linear  waves  by  a  translating  and  oscillating  source-sink  pair  with  zero- 
mean  strength,  (r  =  0.25)  (0<x<ll,0<y<  4.5) 
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Figure  6.15.  Linear  waves  by  a  translating  and  oscillating  source-sink  pair  with  zero- 
mean  strength,  (t  =  0.5)  (n<x<ll,0<y< 4.5) 
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Figure  6.16:  Linear  waves  by  a  translating  and  oscillating  source-sink  pair  with  zero- 
mean  strength,  (t  =  1.0)  (0  <  x  <  11,0  <  y  <  4.5) 


Figure  6.17:  Linear  waves  by  a  translating  and  oscillating  source-sink  pair 
mean  strength,  (t  =  2.0)  (0  <  x  <  11,0  <  y  <  4.5) 


(degree) 


Figure  6.19:  Far  field  linear  waves  by  a  translating  and  oscillating  source-sink  pair 
with  non-zero-mean  strength,  (t  =  1.0) 


Figure  6.20:  Linear  waves  by  a  translating  and  oscillating  source-smk  pair  with  zero- 
mean  strength.  {Fn  =  2  and  t  =  1.0)  (0<x<ll,0<y^  4.5) 


Figure  6.21:  Linear  waves  by  a  translating  and  oscillating  source-sink  pair  with  zero- 
mean  strength.  (F„  =  2  and  t  =  2.0)  (0  <  z  <  11,0  <  y  <  4.5) 
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Figure  6.22:  Constant  phase  lines  of  Kelvin  wake,T  =  0.0 
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Figure  6.23:  Constant  phase  lines,  t  =  0.2 
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Figure  6.24:  Constant  phase  lines,  t  =  0.25_ 
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Figure  6.25:  Constaut  phase  lines,  r  =  0.25+ 


i27 


100.0 

75.0 

^50.0 

25.0 

0.0 


Figure  6.26:  Constant  phase  lines,  t  =  0.5 
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Figure  6.27:  Constant  phase  lines,  t  =  1.0 
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Figure  6.32:  tan“^(|y|/x)  vs.  $  for  third  system  (by  stationary  phase  analysis) 


Figure  6.33:  Bounding-line  angles  qi,  02  and  03  vs  t  (by  stationary  phcise  analysis) 


Figure  6.35:  Linear  waves  due  to  a  translating  and  oscillating  source-sink  pair,  r  = 
0.5,  free  surface  domain:  (0<z<ll,  0<i/<  4.5) 


Figure  6.36:  Nonlinear  waves  due  to  a  translating  ^md  oscillating  source-sink  pair, 
r  =  0.5,  free  surface  domain:  (0<x<ll,  0  <  y  <  4.5) 
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Figure  6.37:  Comparison  of  wave  contour  lines  and  constant  phase  lines  for  linear  | 

waves  due  to  a  translating  and  oscillating  source-sink  pair,  r  =  1.0 
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Figure  6.38:  Unsteady  wake  due  to  a  translating  and  oscillating  source-sink  pair: 
T  =  0.5,  zero-mean  strength,  t  =  12.0 
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0.5,  zero- mean  strength,  <  =  12.5 


Figure  6.40:  Unsteady  wake  due  to  a  translating  and  oscillating  source-sink  pair: 
T  =  0.5,  zero-mean  strength,  t  =  13.0 
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0.5,  zero-mean  strength,  t  =  13.5 


Figure  6.42:  Unsteady  wake  due  to  a  translating  and  oscillating  source-sink  pair: 
r  =  0.6,  zero-mean  strength,  t  =  14.0 
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Figure  6.44:  Wave  elevation  along  inner-angle  wavepacket:  t  =  0.5,  zero-mean 
strength,  12  <  t  <  14) 


Figure  6.45:  Peak  wave  amplitude  within  the  diverging  wave  region  in  nonlinear 
steady  wake 
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Figure  6.46:  Peak  wave  amplitude  within  the  interference  region  in  nonlinear  steady 
wake 
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Figure  6.47:  Peak  wave  amplitude  within  the  diverging  wave  region  in  linear  steady 
wake 


Figure  6.48:  Peak  wave 


amplitude  within  the  interference  region  in  linear  steady 


Figure  6.49:  Peak  wave  amplitude  within  the  inner*angle  wavepacket  in  nonlinear 
imsteady  wake 


Figure  6.50:  Peak  wave  amplitude  within  the  inner-angle  wavepacket  in  linear  un 
steady  wake 
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Figure  6.51:  Peak  wave  amplitude  within  the  inner-angle  wavepacket  in  the  wake  of 
the  cutter  Point  Brower  (experimental  data  of  Brown  et  al.  redrawn  to 
show  the  algebraic  decay  rate) 


CHAPTER  VII 


CONCLUSIONS  AND  SUGGESTIONS  FOR 
FURTHER  RESEARCH 

7.1  Conclusions 

The  desingulaxized  boundary  integral  method  significantly  reduces  the  time  re¬ 
quired  to  compute  the  influence  coeflficient  matrix  of  the  resulting  algebraic  system. 
The  indirect  method  works  better  than  the  direct  method,  especially  for  the  un¬ 
bounded  flows.  The  matrix  is  still  adequately  well  conditioned  for  efficient  iterative 
solutions.  Moreover,  preconditioning  of  the  matrix  greatly  reduces  the  time  to  solve 
the  algebraic  system.  For  the  rases  so  far  investigated,  few  difficulties  of  non-physical 
reflection  from  the  open  computation  boundaries  are  encountered.  Accurate  solu¬ 
tions  can  be  obtained  for  a  large  range  of  desingularization  distances  on  the  order  of 
the  mesh  size.  The  method  is  easy  to  code  and  vectorize  to  achieve  high  performcince 
on  supercomputers.  The  method  has  been  successfully  applied  to  several  nonlinear 
gravity  wave  problems. 

In  the  two-dimensional  upstream  runaway  soliton  problem,  the  fully  nonlinear 
model  predicts  the  runaway  solitons  as  do  the  Boussinesq  amd  fKdV  models.  Our 
numerical  results  agree  with  the  fKdV  results  for  small  disturbances.  However,  for 
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strong  disturbances,  the  fully  nonlinear  model  predicts  larger  solitons  than  the  fKdV. 
These  soliton  are  so  steep  that  they  will  break  sometime  either  upstream  or  down¬ 
stream  of  the  disturbance,  which  have  been  observed  in  experiments.  For  an  equiv¬ 
alent  distribution,  the  free  surface  pressure  generates  significantly  larger  waves  tham 
the  bottom  topography. 

In  the  three-dimensionail  computations  for  the  submerged  spheroid,  the  wave  el¬ 
evation,  the  pressure  distribution  on  the  body  and  the  total  forces  (drag,  lift  and 
moment)  are  calculated.  The  convergence  with  respect  to  the  mesh  size,  desingular- 
ization  distance,  and  energy  conservation  atre  checked.  The  comparison  of  the  drag, 
lift  and  moment  to  other  cadculations  is  reasonably  good  for  Hj D  =  0.245.  For 
HJ D  =  0.16,  our  calculations  predict  surface  piercing  of  the  body.  This  may  be  due 
to  the  lairge  initial  acceleration  or  high  speed  of  the  body  (large  disturbamce).  The 
waves  generated  by  the  spheroid  and  the  relevauit  source-sink  disturbance  are  simi¬ 
lar  if  they  are  sufficiently  submerged.  This  indicates  that  a  simple  disturbance  can 
model  the  body  if  the  far-held  waves  are  the  main  interest.  This  results  in  significant 
computational  savings. 

The  nonlinear  and  linear  computations,  as  well  as  the 'stationary  phase  analysis 
for  the  waves  due  to  the  source-sink  pair  or  surface  pressure  patch,  have  shown  that 
nonlinear  effects  are  not  necessary  in  the  generation  of  the  inner-angle  wavepackets 
while  the  unsteadiness  of  the  flow  can  cause  the  inner-angle  waves.  The  wave  packet 
angle  decreases  with  the  the  increase  of  the  reduced  frequency.  The  inner-angle  wave 
packet  decays  at  the  same  rate  as  the  Kelvin  diverging  waves.  It  appears  that  the 
angle  and  persistence  of  the  inner-angle  wavepacket  behind  the  ship  can  be  explained 
by  the  unsteady  linear  wave  theory.  The  10.8°  inner-angle  observed  in  the  experiment 
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of  the  Coast  Guard  cutter  Point  Brower  corresponds  to  1.0  <  r  <  1.5. 

7.2  Suggestions  for  Further  Research 

Desingularization  greatly  simplifies  the  boundary  integration  and  also  reduces  the 
computation  time.  However,  the  desingularizied  boundary  integral  method  is  still 
an  0{N‘^)  algorithm.  When  preconditioning  of  the  influence  matrix  is  employed,  the 
calculation  of  the  influence  matrix  and  the  velocity  on  the  free  surface  dominates 
the  computation  time  (solving  the  matrix  equation  takes  only  20%  to  30%  of  the 
total  CPU  time)  for  N  «  3000.  The  GMRES  is  also  cin  0{N'^)  algorithm.  The 
recently  developed  Fast  Multipole  Method  for  the  so-called  jV-body  potential  prob¬ 
lem  (Rokhlin  1985,  Greengard  &:  Rokhlin  1988,  Greengard  1990,  and  Korsmeyer 
1990)  seems  promising  to  be  incorporated  with  the  desingularization.  The  multi¬ 
pole  method  evaluates  the  influences  among  the  N  bodies  (mass  particle,  electronic 
charges,  sources,  dipoles,  etc.)  using  multipole  expansions  to  reduce  the  operation 
count  from  0{N^)  to  0{N  log  N).  When  the  multipole  method  is  used  with  an 
iterative  solver  for  a  boundary  value  problem,  it  is  not  necessary  to  form  the  in¬ 
fluence  matrix  explicitly.  This  not  only  reduces  the  memory  storage  requirement 
but  also  reduces  the  operations  for  the  iterative  solver  from  0{N^)  to  0{N  log  N). 
It  is  expected  that  the  computation  time  could  be  greatly  reduced  for  very  large 
scale  calculations  with  the  combination  of  desingularization,  multipole  expansions 
and  preconditioning. 

Parallel  architecture  shows  great  potential  to  speed  up  computations  (Duncan 
1990,  Carriero  &  Gelernter  1989).  There  are  great  potentials  for  parallel  computa¬ 
tions  using  the  desingularized  method.  Parallel  computations  may  become  impera¬ 
tive  if  the  multipole  method  is  used,  since  the  sorting  of  the  particles  in  the  multipole 
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method  requires  many  IF  statements  which  degraule  the  vectorization  of  the  program 
but  not  parallelism. 

Introducing  a  small  number  of  Havelock  sources  into  nonlinear  wave  calculations 
can  improve  the  accuracy  of  the  results  (Scragg  &  Talcott  1990).  The  far-field 
conditions  C2ui  be  better  satisfied  and  a  smzdler  free  surfaoe  computational  domain 
may  be  required.  However,  numerical  instability  in  the  evaluation  of  Havelock  source 
can  sometimes  cause  difficulties.  The  use  of  the  Fourier- Kochin  representation  to 
avoid  these  difficulties  has  been  under  investigation  (Noblesse  &  Hendrix  1990).  For 
unsteady  calculations,  large  memory  storage  (proportional  to  the  simulation  time 
because  of  the  time  convolution  integration)  is  required  to  store  the  influences  due 
to  the  sources  for  each  of  the  previous  time  steps.  Before  this  problem  is  solved, 
large  time  simulations  are  unlikely. 

The  present  work  should  be  extended  to  the  case  of  surface-piercing  bodies  to 
solve  practical  seakeeping  problems,  such  as  the  extreme  motions  of  a  ship  in  a  sea¬ 
way.  This  is  very  important  in  design.  The  effects  of  the  ship  hull  on  the  interaction 
of  bow  and  stern  waves  will  also  allow  a  better  understanding  of  the  inner-angle 
wave  phenomenon.  The  ship  motions  and  the  waves  can  be  solved  together.  How¬ 
ever,  in  addition  to  the  difficulties  associated  with  the  discretization  of  the  wetted 
ship  surface  and  the  free  surface  which  change  with  time,  the  singular  behaviors  of 
the  solution  at  the  intersection  line  of  the  ship  surface  and  the  free  surface  is  the 
most  important  and  difficult  problem.  The  rapid  change  of  the  boundary  conditions 
of  the  two  surfaces  at  the  intersection  line  is  the  cause  of  the  singularity.  In  time¬ 
stepping  computations,  the  specification  of  the  initial  conditions  (or  the  boundary 
conditions  for  the  first  time  step)  has  a  significant  effect  on  the  solution  even  in  the 
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very  simple  wavemaker  problem  (Joo,  Schultz  &  Messiter  1990).  Some  alternative  to 
the  standard  time-stepping  algorithm  used  in  this  study  may  need  to  be  developed. 

Our  conclusion  about  nonlinear  effects  on  the  formation  and  persistence  of  the 
inner-angle  wavepacket  is  based  on  our  limited  nonlinear  wave  calculations  in  which 
the  strength  of  the  disturbance  is  restricted  so  that  no  wave  breaking  occurs.  There¬ 
fore,  the  disturbances  used  may  not  be  sufficiently  strong  to  produce  significantly 
nonlinear  waves  in  the  far  field.  Wave  breaking  was  observed  in  the  cutter  Point 
Brower  experiment.  The  far-field  nonlinearity  may  play  an  important  role,  for  in¬ 
stance  in  making  the  inner-angle  waves  more  isolated  (or  sharper).  Our  method  must 
be  modified  to  continue  computation  after  the  wave  breaking  occurs.  The  main  dif¬ 
ficulty  is  very  similar  to  that  in  the  surface  piercing  body  problem.  The  free  surface 
is  no  longer  smooth  at  the  point  where  the  wave  starts  to  break.  Therefore,  the 
solution  to  the  surface  piercing  body  problem  may  provide  some  insight  in  solving 
the  wave  breaJcing  problem. 
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APPENDIX  A 

Tests  for  Multigrid  Method 


Multigrid  methods  have  been  widely  used  and  proved  effective  in  various  finite 
difference  schemes  (differentiation  type).  Although  the  application  of  multigrid  meth¬ 
ods  to  boundary  integral  methods  (integration  type)  zure  relatively  limited,  it  could 
be  a  good  alternative  to  speed  up  the  computation.  Here,  we  examine  the  application 
of  a  three-level  W  type  multigrid  method  to  our  desingularized  boundary  integral 
method.  Briggs  (1987)  provides  a  good  introduction  to  mutligrid  methods. 

For  an  approximation  x*  to  the  equation  (6.1) 

Ax  =  b.  (A.l) 

we  have  the  residual  vector  r  =  b  — Ax*.  Defining  the  error  as  the  difference  between 
the  exact  solution  to  (6.1)  and  the  approximation,  e  =  x  —  x*,  we  have  the  following 
equation  for  the  error  in  terms  of  the  residual; 


Ae  =  r.  (A.2) 

The  principle  of  multigrid  is  that  diflFerent  frequency  components  of  the  error  are  most 
efficiently  handled  on  different  grids.  Low  frequency  components  are  best  minimized 
on  a  coarse  grid  while  high  frequency  components  are  best  minimized  on  fine  grids. 
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For  ease  of  presentation  we  address  the  case  of  two  grids.  Extension  to  more  grids 
are  readily  made.  In  the  following,  the  superscripts  /  and  c  denote  the  fine  grid 
and  the  coarse  grid  respectively.  Let  x-^  be  the  approximation  to  (6.1)  after  several 
iterations  by  an  iterative  solver  on  the  fine  grid,  then  the  residual  on  the  fine  grid  is 

(A. 3) 

The  residual  on  the  fine  grid  is  interpolated  to  the  coarse  grid  via  a  “restriction” 
operator  R: 

=  Rr^.  (A.4) 

The  error  equation  (A.2)  is  solved  (directly  or  iteratively)  on  the  coarse  grid: 

AV  =  (A.5) 

The  error  e®  is  interpolated  to  the  fine  grid  via  a  “prolongation”  operator  P: 

=  Pe'.  (A.6) 

The  latest  fine  grid  iterate  is  then  corrected  via: 

x^^x^  +  e^.  (A.7) 

The  co2u:se  and  fine  grid  matrices  and  A'  we  set  up  in  an  identical  manner; 
the  only  difference  is  the  truncation  used.  Equation  (6.1)  is  again  solved  on  the 
fine  grid  by  the  iterative  solver  using  the  multigrid-corrected  iterate  as  the  initial 
approximation. 

In  our  computations,  if  the  desired  accuracy  is  not  achieved  within  the  specified  K 
iterations,  the  multigrid  correction  is  applied  and  the  iterative  solver  is  called  again 
for  another  K  iterations.  The  use  of  multigrid  may  affect  the  performance  of  the 
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iterative  solver  used.  The  performance  of  SOR  (successive  over-relaxation),  which 
is  most  widely  used  in  multigrid  methods,  is  not  significantly  affected  by  multigrid. 
However,  GMRES  (generalized  minimal  residual  method — an  extension  of  conjugate 
gradient  method),  which  is  much  more  effective  than  SOR  for  nonsymmetric  matri¬ 
ces,  can  be  degraded  by  multigrid.  Consider  an  approximate  solution  of  the  form 
X*  -i-  Sx.,  where  x*  is  the  initial  guess  and  is  the  correction  to  x*.  Sx  is  a  member 
of  the  Krylov  space  E  =  span(ro,  Afo,  A^r^, ....,  A*“'ro),  in  which  To  =  b  —  Ax*  is 
the  residual  and  fc  (<  A"  N)  is  the  dimension  of  E.  The  GMRES  algorithm  deter¬ 
mines  6x  such  that  the  2-norm  of  the  residual  (i.c.  ||b  —  A(x*  -I-  5x)||  )  is  minimized. 
k  is  also  the  number  of  iterations  conducted,  which  means  that  the  dimension  of  E 
increases  as  the  iteration  number  increases.  R"  is  a  specified  upper  limit  of  k  before 
restarting  of  a  GMRES  cycle.  Usually,  the  performjmce  of  the  GMRES  is  poor  if  K 
is  too  small.  Fig.  A.l  shows  the  the  2-norra  of  the  residual  vs  the  number  of  itera¬ 
tions  for  three  cases  for  a  total  of  100  iterations  without  a  multigrid  correction.  The 
GMRES  is  restarted  every  K  —  b  iterations  in  the  first  case,  iif  =  20  iterations  in 
the  second  case  and  K  =  100  in  the  third  case.  As  czm  be  seen,  after  100  iterations, 
the  third  case  has  a  much  smaller  residual  than  the  other  two  cases.  Also  shown  is 
the  result  for  the  SOR.  The  GMRES  is  better  than  the  SOR  unless  it  is  restarted 
too  frequently  and  therefore  is  degraded. 

Fig.  A.2  shows  the  effect  of  multigrid  on  the  GMRES.  A  coarse  multigrid  correc¬ 
tion  is  applied  every  K  iterations  (no  multigrid  correction  for  K  =  100).  For  K  =  5, 
the  multigrid  reduces  the  residual  faster  in  the  first  few  iterations  (<  20  iterations) 
but  slows  down  the  convergence  significantly  afterward.  For  K  =  20,  the  multigrid 
does  not  help  even  in  the  first  few  iterations  although  the  GMRES  with  =  20 
performs  better  than  with  K  =  5.  Fig.  A.3  shows  the  effects  of  multigrid  on  the 
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SOR.  The  multigrid  accelerates  the  convergence  for  the  SOR.  However,  it  does  not 
help  much  if  a  large  number  of  iterations  are  required  to  achieve  the  desired  accuracy. 
In  our  nonlinear  wave  calculations,  we  usually  require  the  2-norm  of  the  residual  to 
be  less  than  10“®  to  10“^,  which  meams  that  more  than  100  iterations  are  required 
if  the  GMRES  is  used.  Therefore,  the  GMRES  without  multigrid  correction  is  still 
superior. 

Another  reason  why  the  multigrid  correction  may  not  help  in  our  application  of 
the  desingularized  method  is  that  the  integrations  over  the  auxiliary  boundary  must 
be  evaluated  more  carefully  (a  distributed  source  instead  of  isolated  sources  must  be 
used)  and  the  desigularization  distance  should  not  change  for  different  grids  so  that 
the  residual  can  be  restricted  to  the  coarser  grids  and  the  errors  prolonged  up  to  the 
finer  grids.  This  not  only  increases  the  computation  time  in  setting  up  the  matrix 
A  but  may  also  cause  convergence  problems  for  the  desingularized  method  since  the 
ratio  of  the  desigularization  distance  to  the  local  mesh  size  is  no  longer  the  same  on 
the  different  grids. 


Residual:  ||6— /l*a:|| 
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APPENDIX  B 

Preconditioning  of  Equation  (6.1) 


We  seek  another  alternative  to  accelerate  the  convergence:  preconditioning  of 
(6.1).  The  basic  idea  of  preconditioning  is  to  find  an  approximate  matrix  A  to  A  so 
that  A~*  A  is  close  to  the  identity  matrix  or  a  strongly  diagonally  dominated  matrix 
so  that  the  preconditioned  equation: 

A~‘Ax  =  A‘*b  (B.l) 

is  better  conditioned  than  the  original  system  (6.1).  Thus,  it  can  be  solved  by  an 
iterative  solver  with  much  fewer  iterations.  Clearly,  one  requirement  for  A  is  that  it 
be  e^^siiy  invertible  or  decomposed  (z.c.,  a  linear  system  having  A  as  the  coefficient 
matrix  can  be  solved  with  much  less  effort  than  solving  the  original  system).  To 
reduce  the  number  of  iterations,  it  is  also  desirable  that  A  be  close  to  A.  Therefore, 
in  choosing  a  preconditioner,  one  must  select  between  methods  which  usually  perform 
a  large  number  of  “cheap”  iterations  or  a  small  number  of  “expensive”  iterations. 

There  exist  many  preconditioning  methods,  such  as  the  Richardson  method,  the 
incomplete  LU  decomposition  method,  the  modified  LU  decomposition  method,  etc. 
and  the  subroutines  based  on  these  methods  are  available  in  the  NSPCG  (Non- 
symmetric  Preconditioned  Conjugate  Gradient)  package.  Oppe,  Joubert  k  Kincaid 
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(1988)  give  a  brief  background  on  preconditioners.  Usually,  these  methods  are  very 
effective  for  sparse  matrices.  In  general,  for  a  dense,  nonsymmetric  matrix,  it  is  not 
easy  to  find  a  “close”  matrix  or  a  “closer”  r^i.  trix  has  to  be  obtained  at  a  relatively 
higher  expense  in  the  inversion  or  decomposition  of  the  approximate  matrix.  In 
addition,  the  subroutines  in  the  NSPCG  do  the  preconditioning  every  time  when 
they  are  called,  which  may  result  in  unnecessary  overhead.  In  our  nonhnear  wave 
calculations,  matrix  A  is  function  of  time  and  (6.1)  must  be  solved  each  time  step. 
For  large  time  simulations,  a  large  number  of  time  steps  are  necessary.  If  (6.1)  was 
solved  using  some  subroutine  in  the  NSPCG  for  each  time  step,  the  overhead  would 
be  quite  large. 

However,  we  can  easily  find  a  very  close  approximation  A  in  our  problem.  We 
choose  A  as  A  at  f  =  0,  Ag.  Although  the  free  surface  and  the  auxiliau-y  surface 
change  with  time,  Ao  is  still  a  very  good  approximation  to  A.  Therefore  'vc  can 
invert  Ao  once  and  Ao“*  can  be  used  to  precondition  (6.1)  for  every  time  step.  Ag 
can  be  inverted  by  a  direct  method.  Although  the  inversion  of  A,,  is  expensive,  this 
is  done  only  once  for  the  whole  computation  and  the  computation  time  for  inverting 
Ag  is  very  small  compared  to  the  computation  time  saved  in  solving  the  prec  >ndi- 
tioned  systems.  In  GMRES,  most  of  the  com.putation  time  is  used  to  calculate  the 
multiplication  of  the  matrix  4  and  the  vector  x  on  the  left  hand  i/’He  of  (6.1)  in  each 
iteration.  If  preconditioning  is  employed,  it  needs  an  additional  multiplication  of  the 
matrix  and  the  vector  for  the  GMRES  to  calculate  the  left  hand  side  of  (B.l).  This 
overhead  is  also  small  if  the  total  number  of  iterations  can  be  largely  reduced.  The 
preconditioning  matrix  can  also  be  upgraded  after  a  certain  period  of  time. 

The  following  example  shows  how  effective  the  preconditioning  technique  with  the 
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GMRES  is  in  the  nonlinear  wave  computations.  We  compute  the  nonlinear  waves 
generated  by  a  submerged  source-sink  pair  (used  to  model  the  Coast  Guard  cutter) 
starting  to  move  horizontally  from  the  rest.  3367  node  points  on  the  free  surface  and 
3367  isolated  sources  above  the  free  surface  are  used  resulting  a  3367  x  3367  matrix 
A.  The  computation  is  conducted  on  a  Cray  Y-MP.  The  computed  wave  is  shown  in 
Fig.  6.1.  In  this  computation,  1000  time  steps  are  used,  that  is,  we  need  to  solve  (6.1) 
or  (B.l)  for  1000  times.  If  no  preconditioning  is  applied,  for  each  time  step,  GMRES 
takes  an  average  of  120  iterations  (about  18.5  CPU  seconds)  to  converge  in  solving 
(6.1).  So,  the  total  CPU  spent  on  solving  (6.1)  is  about  18500  seconds.  When  the 
preconditioning  is  applied,  GMRES  takes  only  about  5  to  6  iterations  (about  1.75 
seconds)  to  converge  in  solving  (B.l).  The  totad  CPU  for  solving  the  linear  system 
adds  up  to  1900  seconds  including  the  time  for  the  one-time  inversion  of  which 
is  about  150  CPU  seconds.  So,  the  preconditioning  saves  CPU  time  for  the  solution 
of  the  linear  system  by  more  than  80  %. 

Preconditioning  of  the  matrix  also  delays  the  convergence  difficulty  due  to  large 
desingularization  distances  because  the  larger  the  distance,  the  closer  the  matrix  Ao 
is  to  A.  This  allows  the  use  of  a  relatively  large  desingularization  at  the  same  time 
improving  the  accuracy  of  the  numerical  integrations.  The  preconditioning  technique 
can  also  be  used  in  the  iterative  procedure  for  steady  nonlinear  wave  computations 
by  the  method  of  Jensen,  Mi  k  Soding  (1986). 
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